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

3.4.2. Физико-химическая модель

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

рис.3.5

    Расчеты проводились в диапазоне температур 150-400^oС. Интервал более низких температур в работе не рассматривался ввиду его слабого влияния на рудоносность системы. Принятое в модели деление конвективной системы на изотермические ступени показано на рис.3.5. В нисходящей ветви шаг по температуре взят 25^o в интервале 150-200^oС, и 10^o - в интервале 210-400^oС. В восходящей ветви и зоне рудоотложения до температуры 300^oС установлен шаг 5^o, с тем чтобы получить более подробную характеристику происходящего здесь отложения вещества. В результате длина цепочки реакторов достигала в вариантах модели 50 последовательных ступеней. Связь температуры и давления по ступеням на рис.3.5 примерно соответствует течению раствора через приконтактовую зону магматической камеры (см.рис.3.2). Использование в вариантах модели в качестве исходных растворов для восходящей ветви системы составов, полученных на ступенях 350, 380, 400^oС позволяет, в частности, оценить влияние температуры "очага" системы на металлоносность растворов.

    На каждом шаге расчета раствор переносился со ступени на ступень полностью. Скорость взаимодействия с породой задавалась в модели способом, описанным в разделе 2.2 (уравнение 2.9) с параметром подобия (П/В)₁. Величина этого параметра варьировалась в диапазоне от 0,03 до 10 кг/кг, что существенно превышает его природную дисперсию (см. табл. 3.10). Разрастание зоны метасоматизированных пород во времени подчинялось диффузионной модели Д.С.Коржинского (уравнение 2.11). Предельное соотношение для полной переработки базальта (правая часть неравенства 2.12) было принято равным 80 кг породы на кг раствора, что соответствует плотностям базальта и горячей воды 2,9 и =0,7 г/см³ соответственно, и трещинной пористости =5^o/o [Hyndman, Drury, 1976].

Использованные термодинамические данные.

    Для описания термодинамических равновесий в модели использована 15-компонентная мультисистема Н-О-К-Na-Ca-Mg-Fe-Al-Si-C-S-Cl-Cu-Zn-Pb. В модель не включены Mn и Ba, поскольку они не образуют в нисходящей ветви собственных фаз, и входят в метасоматические минералы изоморфно. Для таких соединений нет надежной термодинамической информации, а без этого достоверный расчет поведения Mn и Ba при формировании гидротермальных растворов оказался пока невозможным. Присутствующий в исходном химическом анализе свежего базальта Mn был пересчитан в табл.3.15 совместно с Mg, а Ba - с Ca.

    Мультисистема включала 45 минералов постоянного состава и водный раствор, свойства которого описаны 50 растворенными частицами. Список минералов и частиц дан в табл.3.16. Ввиду того, что породообразующие минералы метасоматитов (актинолит, эпидот, хлорит) представляют собой твердые растворы, в модель включены как конечные, так и промежуточные члены этих растворов, охватывающие типичные для моделируемых объектов интервалы составов [Gillis, Robinson, 1990]. Для промежуточных хлоритов, актинолита и серицита энергии образования были рассчитаны в предположении идеального многопозиционного смешения, для эпидота - по методу [Bird, Helgeson, 1980]. В модели учтены наиболее важные растворенные формы элементов - свободные ионы, хлоридные и гидроксокомплексы, а для макрокомпонентов - также сульфатные и гидрокарбонатные комплексы. Предварительными расчетами [Гричук и др., 1985] было установлено, что в диапазоне моделируемых условий сульфатные, гидросульфидные, гидрокарбонатные комплексы рудных элементов нигде не преобладают, что позволило их не включать в модель. Использованный список комплексных соединений достаточно типичен для задач моделирования гидротермального процесса [Рафальский, 1993].

    Принятым набором минералов (см.табл.3.16) задаются некоторые упрощения физико-химической модели по сравнению с геологической моделью, что необходимо учитывать при интерпретации. Из них наиболее существенны следующие. (а) В мультисистему не включены смектиты из-за недостаточно надежных термодинамических данных. Вследствие этого в низкотемпературной части системы (150-200^oС) возможны завышенная подвижность Al и сдвиг рН в щелочную сторону. (б) Не учтено изоморфное вхождение Zn и Pb в алюмосиликаты. Это имеет следствием завышение мобилизации рудных элементов из метасоматитов в отсутствие собственных сульфидных минералов.

    Термодинамические данные (g_T) для расчета равновесий в мультисистеме взяты из банка термодинамических данных UNITHERM (каф. геохимии МГУ) по состоянию на декабрь 1995 г. Для расчета g_T простых ("базисных") ионов в банке использована модифицированная модель Хелгесона-Киркхама-Флауэрса [Johnson et al., 1992], для комплексных ионов - экстраполяция по уравнению Рыженко. В качестве опорных для экстраполяции свойств комплексных ионов рудных элементов использовались следующие экспериментальные данные: по комплексам Zn - [Плясунов, Иванов, 1990], Cu - [Варьяш, 1991], Pb - [Seward, 1984]. Для H₂S и комплексных соединений Fe были выполнены специальные работы по уточнению их термодинамических характеристик. Опорные значения g_T приведены в Приложении.

    Для расчета коэффициентов активности растворенных частиц использовалось уравнение Дебая-Хюккеля в третьем приближении. Величины параметра были заимствованы из сводки Р.М.Гаррелса и Ч.Л.Крайста [Гаррелс, Крайст, 1968], а для отсутствующих в ней комплексных ионов принимались равными 4,5 [Рафальский, 1993]. Зависимость от температуры коэффициента при линейном члене уравнения взята по [Helgeson, Kirkham, 1974].

Оценка качества термодинамических расчетов.

    Погрешности термодинамических расчетов, связанные с принятыми характеристиками веществ требуют специального обсуждения. Задача определения таких погрешностей - одна из наиболее сложных в современной термодинамике геохимических процессов, и она еще не получила полного и строгого решения. Различные подходы к этой задаче описаны в [Карпов и др., 1976; Дорогокупец, Карпов, 1984].

    Банк UNITHERM не располагает глобальными процедурами согласования данных, применяемыми в банках "ДИАНИК" и THERMINEOS (см., например, [Дорогокупец, Карпов, 1984]). Помещенные в него термодинамические величины должны быть предварительно согласованы по частным системам. Как известно, экспериментальные данные, а, соответственно, и получаемые из них термодинамические константы характеризуются присутствием в них случайных и систематических ошибок. Приводимые в публикациях величины доверительных интервалов в большинстве случаев характеризуют только случайные ошибки экспериментов. Совместная обработка данных многих экспериментальных исследований позволяет выявить и систематические ошибки отдельных серий данных. Если же экспериментальная информация ограничена, возможность присутствия в ней систематических отклонений сохраняется. Влияние их на результаты расчетов в многокомпонентных системах не может быть определено из величин доверительных интервалов.

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

    Для задач, связанных с гидротермальным рудообразованием большую ценность в этом отношении представляет экспериментальная работа Дж.Хемли с соавторами [Hemley et al., 1992], в которой в равновесии с буферными ассоциациями микроклин-мусковит-кварц и пирит-пирротин-магнетит изучена совместная растворимость в хлоридных растворах наиболее важных сульфидов Cu, Zn и Pb. Важным обстоятельством является то, что в работе [Hemley et al., 1992] были достигнуты равновесные состояния систем. Более близкие по химическому составу к предмету настоящего исследования эксперименты в системе "базальт-морская вода" к сожалению, не пригодны для оценки качества расчетов равновесий, поскольку в них не достигаются равновесные состояния даже при длительности опытов около двух лет.

рис.3.6

    С использованием базы данных UNITHERM нами был проведен расчет задачи, имитирующей опыты Дж.Хемли и др. На рис.3.6 показаны результаты экспериментов и расчетов по содержаниям рудных металлов при 300-400^oС и давлении 500 бар. Сравнение их показывает, что расчеты в целом удовлетворительно согласуются с опытами, хотя для всех металлов существует небольшая систематическая ошибка - для Fe и Zn - в сторону завышения, для Cu и Pb - в сторону занижения. В целом по рудным металлам в этой задаче среднеквадратическое отклонение суммарных концентраций составляет 0,26 логарифмические единицы. Эту величину можно использовать в качестве общей приближенной оценки погрешности расчета равновесий для малых компонентов системы. Для макрокомпонентов системы, концентрации которых определяются уравнениями сохранения массы, качество расчета зависит от точности задания валовых составов исходных растворов и пород. Погрешность этих величин не превышает 20 ^o/o относительных.

    Несовершенство этой оценки точности термодинамических расчетов обусловлено в первую очередь тем, что алюмосиликатные минеральные ассоциации в наших модельных расчетах более сложные, чем в опытах Дж.Хемли и др. [Hemley et al., 1992]. К сожалению, надежными опытными данными по составу равновесных растворов в системах, содержащих совместно эпидот, актинолит, хлорит, альбит и сульфиды, а также - другие алюмосиликаты, присутствующие в модельной мультисистеме (см.табл.3.16) мы не располагаем. Однако использованные для проверки опыты проведены в условиях, сходных с гидротермальными растворами океана по температуре, давлению, величинам рН (рН_Т = 4,5) и содержанию хлоридов, и поэтому растворимости сульфидов рудных элементов в них также достаточно близки к наблюдаемым в океане.

    Выполненная проверка не может гарантировать надежность результатов расчетов по индивидуальным формам нахождения элементов, поскольку в силу конкурентного характера комплексообразования возможна взаимная компенсация ошибок в величине суммарной растворимости элемента. Эта компенсация представляет существенную опасность при далекой экстраполяции оценки точности, в особенности - по температуре. Однако опыты Дж.Хемли проведены в том же интервале условий, что и использованные в наших расчетах (300-400^oС, 500 бар против 150-400^oС, Р_{нас.пара} - 500 бар в модельных расчетах). Поэтому оценка точности суммарных концентраций не изменится даже при наличии скомпенсированных ошибок.

Таким образом, общая погрешность расчета равновесий в рассматриваемой системе оценивается нами как 0,1 лог.единиц для концентраций макрокомпонентов системы, и 0,5 лог.единиц для микрокомпонентов.