Геохимия, 2022, T. 67, № 10, стр. 942-960

ВВЕДЕНИЕ

Карбонаты кальция и магния являются наиболее распространенными карбонатами в земной коре. Они представлены соединениями CaCO₃ (тригональный кальцит, ромбический арагонит), CaMg(CO₃)₂ (доломит) и MgCO₃ (магнезит). При высоких давлениях устойчивые на поверхности модификации переходят в высокобарические фазы. Так, тригональный кальцит (cal-I) преобразуется в cal-II с моноклинной решеткой при давлении 1.45 ГПа и, далее, в cal-III при 2.2 ГПа (рис. 1). Приблизительно при давлении 4.5 ГПа и температуре 1273 K доломит разлагается на арагонит + + магнезит (Shirasaka et al., 2002), и магнезит остается основным хранилищем окисленного углерода в мантийных ассоциациях перидотитов и эклогитов (например, Fallon, Green, 1989; Hammouda, 2003; Dasgupta et al., 2004) вплоть до условий высокой температуры (2300 K при давлении 15 ГПа) и давления (до 100 ГПа, Isshiki et al., 2004). Большие объемы карбонатов могут погружаться при субдукции до значительных глубин (например, Bebout, 1995). Ромбический CaCO₃ арагонит как биогенного, так и неорганического происхождения распространен и в условиях поверхности Земли, и до глубин более 20 ГПа (рис. 1).

Рис. 1.

Поля устойчивости различных модификаций CaCO₃ (по Bayarjargal et al., 2018). Римские цифры – обозначение полиморфных модификаций. Цифры в кружках (1) граница арагонит (arg)–кальцит(cal)-II, (2) границы полей устойчивости доломита относительно ассоциации арагонит + магнезит (Shirasaka et al., 2002). Штрих-пунктиром показаны линии геотермического градиента.

Отношения стабильных изотопов в карбонатах, особенно ¹⁸O/¹⁶O и ¹³C/¹²C, определяются температурой кристаллизации, изотопным составом среды образования и часто используются при палеоклиматических реконструкциях (McDermott, 2004; Lea, 2014 и др.). Соотношение “органического” и “карбонатного” углерода отображает баланс углерода в циклическом обмене между океаном, осадками и глубинными породами (Meyer et al., 2013; Wang et al., 2017) и служит для моделирования углеродного круговорота Земли. Отклонения изотопных отношений углерода морских карбонатных отложений связываются с катастрофическими событиями в истории Земли (Kaufman, Knoll, 1995; Hoffman et al., 1998 и др.), в том числе – с периодами глобальных похолоданий и оксигенаций (Meyer et al., 2013). Кроме того, кальцит выступает в качестве референтного материала при калибровках изотопных геотермометров (Chiba et al., 1989), а также является своего рода “полигоном” при оценках влияния давления на изотопное фракционирование (Gillet et al., 1996; Polyakov, Kharlashina, 1994; Polyakov, 1998; Chacko, Deines, 2008). Кальцит дает возможность объединения точных теоретических вычислений для молекулы CO₂ и наиболее представительных измерений изотопного фракционирования кислорода в силикатах (Clayton et al., 1989; Chiba et al., 1989, Clayton, Kieffer, 1991).

Помимо изотопов кислорода и углерода, все большее применение находят “нетрадиционные” изотопные системы карбонатов (⁴⁴Ca/⁴⁰Ca, ²⁶Mg/²⁴Mg). Изотопный состав кальция в морских карбонатах используется при реконструкциях глобальных перемещений Ca (De La Rocha, 2000; Heuser et al., 2005), а изотопный состав Mg – при реконструкциях океанического цикла Mg и континентального выветривания (Tipper et al., 2006; Higgins, Schrag, 2010). Так как потоки Mg и Ca обычно связаны с углеродом, изотопный состав Mg и Ca может отображать и глубинные процессы углеродного обмена (DePaolo, 2004; Huang et al., 2015).

Ввиду широкого применения изотопных систем карбонатов в качестве геохимических трассеров, вплоть до настоящего времени продолжаются исследования по определению и уточнению их факторов изотопного фракционирования (Schmidt et al., 2005; Schauble et al., 2006; Chacko, Deines 2008; Rustad et al., 2010; Schauble et al., 2011; Wang et al., 2017 и др.). При этом все большее внимание уделяется неэмпирическим (ab-initio) расчетам, основанным на определениях β-факторов разных фаз с использованием теории функционала плотности (DFT), которые могут устранить сложности и неопределенности, присущие эмпирическим и/или экспериментальным методам (например, проблемы достижения и доказательства изотопных равновесий).

Цель настоящей работы состоит в определении совокупности β-факторов ¹⁸O/¹⁶O, ¹³C/¹²C, ⁴⁴Ca/⁴⁰Ca, ²⁶Mg/²⁴Mg для Ca и Mg карбонатов (кальцита, магнезита, доломита и арагонита) методом “замороженных фононов” DFT в сравнении с ранее полученными результатами. Применение единого подхода для характеристики разных изотопных систем может нивелировать неопределенности, которые могут обуславливаться применением разных методов (в том числе, разных DFT методов) при оценках факторов изотопного фракционирования. Получение разными методами близких значений β-факторов может свидетельствовать о надежности результатов. Также в работе оценивается влияние давления на изотопное фракционирование ¹⁸O/¹⁶O и ¹³C/¹²C.

РАСЧЕТ β-ФАКТОРОВ

Фракционирование изотопов между фазами A и B определяется величиной α_AB = R_A/R_B, R – отношения содержания изотопов в соответствующих фазах. Используя общепринятые обозначения через величину относительного содержания изотопов, δ, можно записать α_AB = (δ_A + 1000)/(δ_B + 1000). Значение равновесного α_AB можно определить на основе частот колебаний разных изотопологов (и далее рассчитанных приведенных отношений статистических сумм, или “β-факторов”), так что при изотопном равновесии 10³ln α_AB = 1000lnβ_A – – 1000lnβ_B (Bigeleisen, Mayer, 1947). Значения β‑факторов кристаллических фаз в гармоническом приближении:

Частоты ν_q,i вычислены методом “замороженных фононов” (CRYSTAL17, Dovesi et al., 2018) (Приложение) , с использованием полноэлектронных орбиталей гауссового типа и применением метода расширенных ячеек для увеличения Nq и достижения достаточной точности вычисления. Достоверность проведенных вычислений подтверждается сравнением с экспериментальными данными как по кристаллографическим параметрам ячеек (отклонения от параметров решетки, полученных экспериментально, не превышают 1.5%), так и по отклонениям частот (см. ниже).

Частоты колебаний определены для расширенных ячеек, соответствующих симметричным преобразованиям исходных (примитивных) ячеек (Приложение) с увеличением их объема (и количества волновых векторов, учитываемых при суммировании) в 16 раз (для кальцита, магнезита и доломита) и 8 раз (для арагонита), что обеспечивает практическую (с точностью не менее 0.01 при температурах выше 0°С) сходимость определения β-факторов.

Для оценки достоверности расчетов, полученные частоты сопоставлены с экспериментальными данными для изотопологов с природным отношением изотопов (табл. 1).

Спектры всех карбонатов в центре зоны Бриллюэна (q = 0) делятся на две группы (Bottinga, 1968; Deines, 2004; Valenzano et al., 2007). Частоты больше ≈700 см^–1 определяют “внутренние” колебательные моды CO₃, а частоты меньше ≈400 см^–1 соответствуют “внешним”, относительно группы CO₃, модам. Максимальные расхождения между вычисленными и экспериментальными значениями частот ранее отмечалось особенно для “внешних” низкочастотных колебаний (например, для кальцита, отклонение низкочастотных колебательных мод A_2u на 33.7 cm^–1 и E_u на 22.4 cm^–1), что связывается или с проблемами экспериментальных определений при низких длинах волн, или с влиянием дисперсионных взаимодействий (Valenzano et al., 2007). Максимальные расхождения колебательных мод арагонита (табл. 1Б) объясняется их нестабильностью (“мягкие” моды, Carteret et al., 2013) и выраженной зависимостью данных мод от объема решетки. Для удобства сопоставлений с ранее полученными результатами приведено среднеквадратичное отклонение вычисленных частот (Δ_MSWD) относительно экспериментальных данных (табл. 1). Без учета аномальных значений Δ_MSWD составляет 5.3 см^–1 (кальцит), 6.9 см^–1 (магнезит) и 9.7 см^–1 (доломит). Среднеквадратичное отклонение вычисленных частот арагонита составляет 7.6 см^–1 (11.5 см^–1 с учетом аномальных мод). Более показательно сопоставление с учетом регрессионной зависимости между вычисленными и экспериментальными (по ИК- или Рамановским спектрам) частотами $\omega _{i}^{{\exp }}$ = SFω_i, где SF – так называемый масштабный фактор – коэффициент линейной регрессии между вычисленными и наблюдаемыми частотами (Schauble, Young, 2021). Масштабный фактор определяет систематическое отклонение вычисленных значений частот относительно экспериментальных и в ряде случаев используется в качестве поправки вычисленных частот и, соответственно, величин β-факторов (см. ниже). Полученные SF отклоняются от единицы не более чем на 0.1–0.4% (табл. 1), а среднеквадратичное отклонение отношений ${{{{\omega }_{i}}} \mathord{\left/ {\vphantom {{{{\omega }_{i}}} {\omega _{i}^{{\exp }}}}} \right. \kern-0em} {\omega _{i}^{{\exp }}}}$ от единицы составляет 0.015–0.05.

Для определения β-факторов и их температурных зависимостей 1) для каждой фазы вычислены изменения частот при изотопных замещениях, 2) по формуле (1) для температур от 0 до 1500°С с шагом 10°C определены 1000 ln β, которые интерполированы с использованием полинома 1000 ln β = =  ax + bx²+ cx³, x = 10^6/T² (табл. 2). Ошибка аппроксимации кубическим полиномом во всех случаях не превышает 0.02%о при температурах выше 100°C и 0.03‰ при температурах выше 0°C (R² > 0.999). По представленным данным можно определить факторы изотопного фракционирования между карбонатами, карбонатами и H₂O–CO₂ флюидами, а также другими фазами с уже известными β-факторами.

Влияние отклонения вычисленных значений частот от экспериментальных спектров на β-факторы можно оценить с помощью масштабного фактора SF. Коэффициенты a, b и c кубического полинома 1000lnβ пропорциональны начальным моментам функции распределения по фононным частотам 2i (Поляков и др., 2019; Shiryaev et al., 2020), поэтому при пропорциональном изменении всех частот на фактор SF, a^SF = aSF², b^SF = = bSF⁴, c^SF = cSF⁶. Полученные в работе значения SF кальцита, магнезита, доломита (табл. 1А) составляют 1.0005–1.0020 и не приводят к сколько-нибудь значительному изменению величин β (табл. 2). Скорректированные с учетом масштабного фактора значения 1000lnβ по кислороду кальцита отличаются на 0.4 при Т = 0°С и на 0.08 при 500°С. Для углерода соответствующие изменения составляют 0.7 при Т = 0°С и 0.15 при 500°С. Для арагонита (табл. 1Б) SF = 0.996 и соответствующие отличия составляют 0.7–0.15 (0–500°C, 1000 ln β¹⁸O) и 1.4–0.25 (0–500°C, 1000 ln β¹³C).

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Сравнение с расчетами β-факторов, выполненными полуэмпирическими методами и “из первых принципов” и изотопным фракционированием в экспериментах

В табл. 2 для сравнения приведены данные расчетов β-факторов карбонатов с использованием значений экспериментально наблюдаемых частот (согласование силовых постоянных с экспериментальными частотами: Chacko et al., 1991; Deines, 2004; Chacko, Deines, 2008, применение унифицированных правил при оценках изотопных сдвигов частот: Clayton, Kieffer, 1991 – далее условно для краткости “полуэмпирические” методы определения β-факторов), а также неэмпирическими (“ab-initio”) методами с применением базиса плоских волн (DFPT – теории возмущений функционала плотности, Schauble et al., 2006; Schauble, 2011) и других DFT методов (Rustad et al., 2010). Учитывая вероятное взаимодействие с углекислотой при образовании карбонатов, в таблице также приведены значения $1000{\text{ln}}{{\beta }_{{{\text{C}}{{{\text{O}}}_{{\text{2}}}}}}}.$

Полученные результаты во многих случаях близки ab-initio вычислениям DFPT (Schauble et al., 2006). Максимальные расхождения проявляются для магнезита (1000 ln β¹⁸O¹¹ и 1000 ln β¹³C при температурах 0–200°С отличаются, соответственно на 0.3–0.2 и 1.9–1.0, сплошные линии на рис. 2). Для кальцита результаты отличаются на 0.8–0.2 (0–200°C) (1000 ln β¹⁸O) и на 0.5 (0°C) – 0.2 (200°C) (1000 ln β¹³C); для доломита в этом же интервале температуры – на 1.1–0.5 (1000 ln β¹⁸O) и 0.1–0.2 (1000 ln β¹³C), для арагонита – на 0.5–0.1 (1000 ln β¹⁸O) и 0.9–0.4 (1000 ln β¹³C). Применение единого масштабного множителя (Schauble, Young, 2021) дает несколько большие значения 1000 ln β¹⁸O (отличие в интервале 0–200°C составляет 0.9–0.6). Результаты полуэмпирических расчетов в том же интервале отклоняются (пунктирные линии на рис. 2) для магнезита на 1–2, (и до более высоких температур), доломита 0.10–0.08 и арагонита 1.1–0.7. 1000 ln β¹⁸O кальцита при температурах до 100°C отличаются более существенно (до нескольких единиц). Все β-факторы по углероду, определенные полуэмпирическим методом (Deines, 2004), существенно отличаются от результатов расчетов из первых принципов (штриховые линии на рис. 2). Расхождение может обуславливаться выбором частот (ν₃, Mironenko et al., 2018).

Рис. 2.

Сравнение 1000 ln β карбонатов, полученных DFT и полуэмпирическими методами. Δ¹⁸O и Δ¹³С по вертикальным осям – отличия 1000 ln β¹⁸O и 1000 ln β¹³C результатов настоящей работы и другими данными: сплошная линия на всех графиках (Schauble et al., 2006); SY (Schauble, Young, 2021); CD (Chacko, Deines, 2008); CM (Chacko et al., 1991); CK (Clayton, Kieffer, 1991); Dei (Deines, 2004). Данные настоящей работы приведены с учетом масштабного множителя. Области совпадения результатов (Δ ≈ 0) затемнены.

Согласование расчетов DFT с экспериментальными данными может быть критичным при оценке достоверности различных методов. Возможные ошибки и расхождения методов DFT обуславливаются неизбежными приближениями при выборе обменно-корреляционных функционалов взаимодействия электронов, базисных функций для описания атомов, численными ошибками при расчетах. Ошибки экспериментальных результатов нередко вызываются проблемами достижения и контроля изотопного равновесия, возможным перераспределением изотопов при охлаждении, влиянием других физико-химических условий экспериментов (давление, составы реагентов). На рис. 3 показаны результаты сравнения факторов изотопного фракционирования кислорода и углерода в системе карбонат-CO₂, полученными методами DFT, расчетами по модельным спектрам, с наиболее представительными результатами экспериментов прямого обмена, а также осаждения карбонатов (кальцита). Расчетные значения факторов фракционирования ¹⁸O/¹⁶O между CO₂ и кальцитом практически не отличается от “эталонных” экспериментальных результатов (Chacko et al., 1991). Результаты других экспериментов отклоняются более существенно (более чем на +1.5‰ при T 350–610°C, O’Neil, Epstein, 1966, –1…–1.5‰, T 500–900°C, Scheele, Hoefs, 1992). Можно отметить, что результаты вычислений по модельным спектрам (Поляков, 2008) также соответствуют результатам DFT. Известные экспериментальные определения фракционирования ¹⁸O/¹⁶O между CO₂ и доломитом недостаточно представительны. При T 350–610°C отклонения составляют до +1‰ (O’Neil, Epstein, 1966), при более низких температурах – около ±0.5‰ (Horita, 2014). Известные экспериментальные данные по обмену ¹³C/¹²C между карбонатами и CO₂ в целом удовлетворительно согласуются с расчетными (рис. 3в–3д). Выраженные отличия экспериментов (до ±1‰) наблюдаются при повышенных температурах (>500°C) при обмене кальцит-CO₂ (Scheele, Hoefs, 1992).

Рис. 3.

Фракционирование изотопов ¹⁸O/¹⁶O (а, б) и ¹³C/¹²C (в, г, д) между карбонатами и CO₂.

Теоретические факторы фракционирования (DFT – сплошные, на основе модельных спектров – пунктирные линии) рассчитаны относительно $1000{\text{ln}}{{\beta }_{{{\text{C}}{{{\text{O}}}_{2}}}}}$ (Chacko et al., 1991). Kr – настоящая работа, DFT – усредненные результаты, SGE (Schauble et al., 2006), остальные обозначения – см. рис. 2. Экспериментальные данные (точечные линии): Hor (Horita, 2014); NeEp (O’Neil, Epstein, 1966), ScHo (Scheele, Hoefs 1992), RGM (Romanek et al., 1992), Ros (Rosenbaum, 1994, залитый треугольник), JL (Jiménez-López et al., 2001, залитый кружок), Pol (Поляков, 2008). PMA – (база данных GEOCHEQ_ISOTOPE, Поляков и др., 2021). Тройная сплошная линия (б) – совпадающие результаты CO₂-доломит (Kr, DFT, CD). Затемненные квадраты (б) – Hor (показаны результаты при 100 и 150°С, данные для 250°С существенно варьируют, вероятно их искажение, Horita, 2014). В интервале 20–60°С экспериментальные определения $1000{\text{ln}}{{\alpha }_{{{\text{C}}{{{\text{O}}}_{{\text{2}}}}}}}$ при осаждении кальцита (Emrich et al., 1970) практически совпадают с результатами DFT расчетов (на рис. 3в линии сливаются).

Сравнение вычисленного фракционирования ¹⁸O/¹⁶O между карбонатами (кальцитом, доломитом) и H₂O с экспериментальными данными представлено на рис. 4. При температурах меньше 100°С вычисленные 1000 ln  α¹⁸O кальцит-H₂O (рис. 4а) практически совпадают с экспериментальными результатами (O’Neil et al., 1969), а при температурах менее 20°С – с результатами (Horita, Clayton, 2007). Значения 1000ln α¹⁸O доломит-H₂O отличаются незначительно (0.1–0.2‰) от оценок, основанных на комбинации экспериментальных результатов (Horita, 2014). При T 200–350°C (рис. 4б) отклонение оценок с применением β-факторов жидкой воды (Schauble, Young, 2021) дает отклонения около –0.5‰ относительно экспериментальных определений (O’Neil et al., 1969) и при T > 350°C расхождения оценок фактора фракционирования увеличиваются. При расчетах с использованием β-факторов водяного пара отклонения уменьшаются практически до нуля (приближение к критическому состоянию H₂O). Аналогично, при T до 300–350°С вычисленные значения 1000 ln α¹⁸O приближается к экспериментальным при равновесии доломита с водяным паром, а при T > 350°C – при равновесии доломита с H₂O в надкритическом состоянии (рис. 4б).

Рис. 4.

Фракционирование изотопов ¹⁸O/¹⁶O (а: 20–100°C; б: 200–600°C) между карбонатами и H₂O.

Теоретические факторы фракционирования (кальцит-H₂O – сплошные, доломит-H₂O – штрихпунктирные линии) рассчитаны относительно $1000{\text{ln}}{{\beta }_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}}$ (Schauble, Young, 2021: жидкой воды – H₂OL, водяного пара – H₂OV, надкритической фазы – H₂Osc).

Экспериментальные данные по кальциту (пунктирные линии): NCM (O’Neil et al., 1969), KNH (Kim et al., 2007), HC (Horita, Clayton, 2007) и доломиту (штриховая линия): Hor (Horita, 2014), NC (Norton, Clayton, 1966). Результаты Kim, O’Neil, 1997 совпадают с HC и на рисунке не показаны.

В целом можно отметить, что расхождение оценок, полученных “из первых принципов” в ряде случаев меньше, чем отклонения различных экспериментальных данных между собой. Минимальное расхождение получено по фракционированию ¹⁸O/¹⁶O в системе кальцит-CO₂ (эксперименты Chacko et al., 1991), более выражены отличия факторов фракционирования ¹³C/¹²C. На определение величины изотопного фракционирования карбонатов с H₂O существенно влияет различие β-факторов в зависимости от условий (жидкость, пар, надкритическое состояние).

Вычисленные β-факторы “нетрадиционных” изотопов магния (1000 ln β²⁶Mg) для магнезита практически совпадают с результатами DFPT (Schauble, 2011), для доломита отличия составляют 0.5–0.2 (0–200°C). Отличие результатов других известных DFT вычислений (Rustad et al., 2010) составляет для магнезита 1.4–0.5 (1000 ln β¹⁸O) и 3.1–1.1 (1000lnβ¹³C) для доломита. Вычисленные в том же интервале температур 1000lnβ⁴⁴Ca отличаются на 2.5–0.9 (кальцит, Rustad et al., 2010) и 0.8–0.2 (доломит, там же).

Фракционирование магнезит – кальцит определяется выражениями (по данным табл. 2, с учетом масштабного множителя):

(2)

$\begin{gathered} {{10}^{3}}\ln {{\alpha }^{{18}}}{{{\text{O}}}_{{mgs}}}_{{{\text{ - }}cal}} = 1.01010x - \\ - \,\,0.02787{{x}^{2}} + 0.00078{{x}^{3}}, \\ {{10}^{3}}\ln {{\alpha }^{{13}}}{{{\text{C}}}_{{mgs}}}_{{{\text{ - }}cal}} = {\text{ }}1.14764x - \\ - \,\,0.08387{{x}^{2}} + 0.00294{{x}^{3}}. \\ \end{gathered} $

Согласно (2), при 0°С равновесное изотопное фракционирование ¹⁸O/¹⁶O между магнезитом и кальцитом достигает 10.4‰ (4.0‰ при 200°С), фракционирование ¹³С/¹²С при 0°С составляет 7.4‰ (3.7‰ при 200°С). Для сравнения, по данным расчетов методом DFPT с применением базисных функций в виде плоских волн и псевдопотенциалов (Schauble et al., 2006) 10³ln α¹⁸O_mgs-cal(0°С) = = 11.4‰ , 10³lnα¹³C_mgs-cal(0°C) = 6.0%o. Таким образом отличия результатов, полученных различными методами DFT, при низких температурах составляет 1‰ (¹⁸O/¹⁶O) и 1.4‰ (¹³C/¹²C). Результаты, полученные с использованием экспериментальных частот для определения силовых постоянных (Chacko, Deines, 2008; Deines, 2004) отличаются от результатов DFT более существенно: 10³ln α¹⁸O_mgs-cal(0°С) = 14.0‰, 10³ ln α¹³C_mgs-cal(0°С) = = 7.1‰. Отличия температур, вычисленных с применением разных калибровок в диапазоне 0–200°С, возрастают и особенно существенны при сравнениях результатов ab-initio расчетов с полуэмпирическими калибровками (Deines, 2004; Chacko, Deines, 2008), что обуславливается неоднозначным выбором частот при моделировании спектров.

При изоморфном вхождении Mg в решетку кальцита факторы изотопного фракционирования магнезит-кальцит в целом уменьшаются. Для корректировки β-факторов кальцита можно воспользоваться результатами определения 1000lnβ в зависимости от содержания Mg (рис. 5):

(3)

$\begin{gathered} 1000\ln {{\beta }_{{{\text{cal}}}}}\left( {x,{\text{Mg}}} \right) = 1000\ln {{\beta }_{{{\text{cal}}}}}\left( x \right) + \\ + \,\,\Delta 1000\ln {{\beta }_{{{\text{cal}}}}}\left( {{\text{Mg}}} \right) = \left( {a + \Delta a} \right)x + \left( {b + \Delta b} \right){{x}^{2}} + \\ {\text{ + }}\,\,\left( {c + \Delta c} \right){{x}^{3}}, \\ \end{gathered} $

где a, b, c – коэффициенты полиномиальной зависимости 1000lnβ без учета вхождения Mg в решетку кальцита (табл. 2), Δa, Δb, Δс – изменение указанных коэффициентов за счет вхождения Mg (рис. 3). Необходимо отметить, что отклонение от идеальности изоморфизма приводит к существенному отклонению от линейного соотношения 1000lnβ vs Mg_cal.

Рис. 5.

1000lnβ кальцита в зависимости от содержания Mg. Δa, Δb, Δc – изменение соответствующих коэффициентов полиномов температурной зависимости 1000 ln β кальцита (табл. 2) при мольной доле Mg в кальците, m.

Фракционирование изотопов кислорода и углерода между доломитом и кальцитом. Изотопный термометр доломит-кальцит – один из наиболее применяемых при изучении происхождения карбонатных пород. Термометр первоначально был откалиброван эмпирически (Goldsmith, Graph, 1958) на основе содержания Mg в кальците и кальцит-доломитового сольвуса. Экспериментальная калибровка (Sheppard, Schwartz, 1970) значительно отличается от эмпирической вследствие возможного перераспределения изотопов при охлаждении. Фракционирование изотопов в случае, если кальцит представлен чистым CaCO₃, по результатам проведенных в настоящей работе расчетов дается выражениями:

(4)

$\begin{gathered} {{10}^{3}}\ln {{\alpha }^{{18}}}{{O}_{{dol}}}_{{{\text{ - }}cal}} = 0.51533x - 0.01634{{x}^{2}} + \\ + \,\,0.00043{{x}^{3}}, \\ {{10}^{3}}\ln {{\alpha }^{{13}}}{{C}_{{dol}}}_{{{\text{ - }}cal}} = 0.43095x - 0.03268{{x}^{2}} + \\ + {\text{ }}0.00100{{x}^{3}}. \\ \end{gathered} $

При этом 10³ln α¹⁸O_dol-cal (25°С) = 4.3‰, 10³ln α¹³C_dol-cal (25°С) = 2.1‰, что очень близко результатам других определений (Schauble et al., 2006 при 0°С отличаются на 0.3 и на 0.2‰ при 200°С), так что вычисленные в диапазоне 0–200°С температуры отличаются не более чем на 10–30°C (рис. 6).

Рис. 6.

Отклонение температуры по доломит-кальцитовому геотермометру относительно результатов настоящей работы. Одинарные линии – по распределению изотопов кислорода, двойные линии – по распределению изотопов углерода. Обозначения – см. рис. 2. Пунктиром ограничена область ошибки при измерении величин δ¹⁸O с точностью ±0.05‰. Данные настоящей работы приведены с учетом масштабного множителя.

Изоморфное замещение Ca на Mg в кальците, наблюдаемое как в условиях накопления осадков, так и при повышенных температурах, существенно влияет на изотопное фракционирование и должно учитываться при расчетах (Jiménez-López et al., 2004; Chacko, Deines, 2008). Если твердый раствор CaCO₃–MgCO₃ идеальный, то справедливо соотношение:

(5)

$\begin{gathered} {{10}^{3}}\ln \alpha \left( {dol{\text{ - }}cal} \right) = 1000\ln {{\beta }_{{dol}}} - 1000\ln {{\beta }_{{cal}}} - \\ - \,\,\left( {1000\ln {{\beta }_{{mgs}}} - 1000\ln {{\beta }_{{cal}}}} \right)m, \\ \end{gathered} $

где m – мольная доля магния в кальците. Результаты близки экспериментальным значениям изотопного фракционирования при синтезе Mg-кальцита (при 25°C 10³lnα¹⁸O возрастает 0.17‰ на моль MgCO₃, Jiménez-López et al., 2004) и составляют при той же температуре 0.13‰ на моль MgC-O₃. Учитывая, особенно с увеличением содержания Mg, отклонение от идеальности, при вычислениях факторов изотопного фракционирования (рис. 5) вместо линейного соотношения (5) справедливы более сложные зависимости:

10³lnα(dol-cal) = 1000 ln β_dol – 1000 ln β_cal(x,Mg), где 1000lnβ_cal(x,Mg) определяется соотношением (3).

Согласно полученным в настоящей работе данным, равновесное фракционирование изотопов магния между магнезитом и доломитом при низкой T (0°С) составляет –0.98‰, а при 200°C –0.32‰. Другие оценки составляют, соответственно –1.36, –0.46‰ (Schauble et al., 2006) и ‒2.64, –0.93‰ (Rustad et al., 2010). Фракционирование изотопов кальция по нашим результатам составляет 0.42‰ (cal-dol), 4.55‰ (cal-arg) при 0°С и, соответственно 0.13 и 1.60 при 200°С. Другие оценки составляют 1.27 и 0.44 (cal-dol 0 и 200°С, Wang et al., 2017), 2.42, 0.85 (Rustad et al., 2010). В паре cal-arg разброс менее значительный (4.55 vs 3.33 Wang et al., 2017 при 0°С).

Влияние давления на β-факторы карбонатов

Влияние давления на β-факторы определяется посредством равенства:

(6)

${{\left( {\frac{{\partial \ln \beta }}{{\partial P}}} \right)}_{T}} = {{\left( {\frac{{\partial \ln \beta }}{{\partial V}}} \right)}_{T}}{{\left( {\frac{{\partial V}}{{\partial P}}} \right)}_{T}} \approx - \frac{1}{{{{K}_{T}}}}{{\left( {\frac{{\Delta \ln \beta }}{{\Delta \ln V}}} \right)}_{T}},$

где K_T – изотермический модуль объемной упругости, V – объем ячейки, P – давление. Величины β-фактора определены при разных объемах кристаллической решетки и выражены в виде поправок к коэффициентам полиномиального разложения 1000lnβ (табл. 3), что согласуется с расчетами влияния P по изменениям частот колебаний (Polyakov, Kharlashina, 1994; Polyakov, 1998). В качестве оценок K_T приняты значения, рассчитанные в рамках квазигармонического приближения (CRYSTAL17, Erba, 2014).

Таблица 3.  

Влияние давления на β-факторы (квазигармоническое приближение)

Минерал	Модуль упругости		Δ1000 ln β/ΔP
	KT, ГПа	KT′	aP	bP	cP
18O/16O
Кальцит	75.40	4.20	0.126660	–0.00557	0.000171
	Polyakov, Kharlashina, 1994*		0.1004	–0.0032	–
	Поляков, 2008**		0.14700	–0.00624	0.00000
Магнезит	105.98	4.59	0.099316	–0.00362	0.000101
Доломит	100.25	4.90	0.105119	–0.0046	0.000144
Арагонит	71.87	5.03	0.117319	–0.00425	0.000128
13C/12C
Кальцит	–	–	0.136937	–0.01129	0.000419
	Поляков, 2008**		0.26835	–0.01256	–
Магнезит	–	–	0.103667	–0.00743	0.000248
Доломит	–	–	0.116308	–0.01007	0.000388
Арагонит	–	–	0.127052	–0.01122	0.000462

Примечания. Поправка 1000 ln β на давление Δ1000 ln β ≈ P(ГПа)*(a_P·x + b_P·x² + c_P·x³), где x = 10⁶/T² (K^–2). Параметры уравнения состояния (EoS) карбонатов в QHA приближении – см. табл. 4.  * Вычисления с использованием параметров Грюнайзена для EoS 3-го порядка. ** Пересчитаны к кубическому полиному (P = 0–3 ГПа).

Таблица 4.  

Уравнения состояния (EoS) карбонатов

Фаза	V0, Å3	E, a.u.	KT, ГПа	KT′	Метод	Источник
Кальцит	121.83	–1883.1557	75.40	4.20	DFT	Наст. работа
	122.60	–	73.46	4.00*	Эксп.	ReAn99
Магнезит	91.59	–928.1875	105.98	4.59	DFT	Наст. работа
	93.18	–	97.1	5.44	Эксп.	LFOKF08
	93.18	–	108.4;110.3;107.2	4.00*	Эксп.	LFOKF08
	93.07	–	108	5.00	Эксп.	FGK02
Доломит	105.20	–1405.6751	100.25	4.90	DFT	Наст. работа
	106.74	–	94	4*	Эксп.	RoRe92
Арагонит	228.63	–3766.2839	71.87	5.03	DFT	Наст. работа
	226.71	–	65.24	4.95	Эксп.	LZCWQ15
	232.50	–	66.09	4.64	DFT	HLHXL17
	227.11	–	65.78	5.10	Эксп.	LSGBD17
	227.11	–	67.08	4.74	Эксп.	LSGBD17

Примечания. Вычисления по алгоритму Erba (2014). * Экспериментальные данные согласованы со значением $K_{0}^{'}$ = 4. LFOKF08 (Litasov et al., 2008), RoRe92 (Ross, Reeder, 1992), ReAn99 (Redfern, Angel 1999), FGK02 (Fiquet et al., 2002), LZCWQ15 (Li et al., 2015), HLHXL17 (Huang et al., 2017), LSGBD17 (Litasov et al., 2017).

Увеличение величины β-фактора в результате компрессии может компенсироваться его уменьшением при тепловом расширении:

(7)

${{\left( {\frac{{\partial \ln \beta }}{{\partial T}}} \right)}_{P}} = {{\left( {\frac{{\partial \ln \beta }}{{\partial V}}} \right)}_{P}}{{\left( {\frac{{\partial V}}{{\partial T}}} \right)}_{P}} \approx {{\alpha }_{T}}{{\left( {\frac{{\Delta \ln \beta }}{{\Delta \ln V}}} \right)}_{P}},$

где α_T – коэффициент объемного теплового расширения (Δβ = 0 если α_TΔT = ΔP/K_T).

Отклонения от линейности зависимости β(V)_T пренебрежимо малы, так что ∂lnβ/∂V постоянные (для каждого карбоната) величины (пример β(V) кальцита при разных температурах показан на рис. 7), K_T меняется в зависимости от давления (учтена первая производная по давлению: K_T ≈ K_{T, 0} + $K_{{T,0}}^{'}P$), α_T зависит от температуры (α_T ≈ α_T + α_TT).

Рис. 7.

Изменение β-фактора в зависимости от объема решетки кальцита (V).

При давлениях до 1.5 ГПа, ограничивающих поле стабильности кальцита-I, максимальный вычисленный эффект давления составляет 1.4‰ по кислороду (1.5‰ по углероду). При более высоких давлениях кальцит-I становится неустойчивым и переходит в высокобарическую фазу кальцит-II, что сопровождается увеличением β (табл. 2). Возможный геотермический градиент ограничивает изменение β-факторов под влиянием давления. При давлении 2 ГПа, изменение 1000 ln β достигает 2.0, 0.8 и 0.3 (T = 0, 200, 500°C), но учитывая возможные геотермические градиенты, может превышать 0.5–1.0 только в условиях сверхдавлений (при геотермическом градиенте менее ≈10°/Км). Выражения зависимости 1000 ln β от давления для Ca–Mg карбонатов представлены в табл. 3.

ЗАКЛЮЧЕНИЕ

Методом “замороженных фононов” теории функционала плотности (DFT) с использованием полноэлектронных базисов впервые определены совокупности β-факторов (¹⁸O/¹⁶O, ¹³С/¹²С, ²⁶Mg/²⁴Mg и ⁴⁴Ca/⁴⁰Ca) Ca-Mg карбонатов в гармоническом и квазигармоническом приближениях при температурах от 0 до 1500°C. Для достижения представительности волновых векторов при суммировании применен метод расширенных ячеек (с увеличение объема в 16 раз для кальцита, магнезита и доломита и в 8 раз для арагонита).

Результаты определения факторов изотопного фракционирования, полученные различными методами DFT, отличаются вследствие применения разных функций для описания атомов (наборы базисов в виде плоских волн или линейные комбинации атомных орбиталей в виде гауссианов, полноэлектронные или с псевдопотенциалами), различных функционалов для описания взаимодействия атомов и других параметров. Наши оценки во многом согласуются с определениями β-факторов методом теории возмущений функционала плотности (DFPT, Schauble et al., 2006). Например, кислородные β-факторы кальцита и магнезита при 25°С отличаются, соответственно на 0.6 и 0.3, а при 300°С – на 0.1 и 0.2. Различие углеродных β-факторов составляет для кальцита 0.5 (25°С) и 0.2 (300°С); для магнезита, соответственно 1.7 и 0.7. Существенные расхождения наблюдаются с некоторыми другими оценками β-факторов, особенно по углероду (Deines, 2004): при 25°С расхождения доходят до 6.3 (кальцит), 6.6, 6.4 (магнезит, доломит), 4.7 (арагонит). Применение “универсального” масштабного фактора (SF = 1.043, Schauble, Young, 2021) приводит к увеличению 1000 ln β¹⁸O на 0.9–1.0 при 0°C (0.5–0.6 при 200°C).

Автор признателен рецензентам Е.О. Дубининой и В.Б. Полякову за конструктивные замечания, которые способствовали дополнению и исправлению работы.

Работа выполнена при поддержке гранта Российского научного фонда № 22-27-00275, https:// rscf.ru/project/22-27-00275/

ПРИЛОЖЕНИЕ

Частоты колебаний вычислены с использованием метода “замороженных фононов” (CRYSTAL17, Dovesi et al., 2018) с применением набора гауссовых полноэлектронных базисов и гибридного функционала B3LYP после оптимизации геометрических параметров решетки, включая размеры и положения всех атомов. В методе “замороженных фононов” энергия фонона вычисляется как функция амплитуд смещений через различия энергий решетки со смещениями и идеальной решетки. Нахождение частот фононов в точке Г(0,0,0) сводится к диагонализации взвешенной на массы атомов матрицы вторых производных энергии по смещениям атомов:

Применялись наборы базисов, Mg 8-511d1G, Ca 86-511d21G, C 6-311d11G, O 8-411d11G(1), (http://www.crystal.unito.it). Точность вычислений контролировалась параметром TOLINTEG 8 8 8 9 27, критерий сходимости энергии самосогласованного поля (SCF) установлен на уровне 10^–12 Хартри. Узлы суммирования обратной решетки (6, 6) соответствует 32 (кальцит, магнезит), 40 (доломит) и 112 (арагонит) независимым k-векторам неприводимой части зоны Бриллюэна. Процесс оптимизации геометрии решетки включал релаксацию координат ядер атомов и параметров решетки по квази-ньютоновскому алгоритму.

Значения вычисленных частот колебаний зависят от выбора q-векторов, по которым проводится суммирование в пределах зоны Бриллюэна. Точность средних значений в выражении (1) улучшается в том числе, при расчетах по расширенным ячейкам (“суперячейкам”), которые получаются путем целочисленных линейных преобразований исходных единичных ячеек. Волновой вектор q расширенной ячейки является вектором обратной решетки, линейные размеры которой должны быть не менее чем 2π/|q|. Определители матриц преобразований det{L} соответствуют количеству волновых векторов при суммировании в (1), т.е. Nq. Для оценки влияний Nq на результаты (β-факторы) для карбонатов частоты колебаний вычислены по расширенным ячейкам разных объемов, соответствующим Nq = 1, 2, 4, 8, 16 и 27. Сходимость (с точностью <0.01‰) достигнута при Nq = 16 (Nq = 8 для арагонита). Соответствующие матрицы преобразования, сохраняющие симметрии исходных ячеек {100/010/001} (единичная ячейка), {211/121/112}(Nq = 4), {200/020/002} (Nq = = 8), {0-2-2/-20-2/-2-20} (Nq = 16), and {300/030/003} (Nq = 27). Такие расширенные ячейки при суммировании эквивалентны 1, 2, 4, 8, 16, и 27 k-точкам первой зоны Бриллюэна (Γ-точке) по схеме Монкхерста–Пэка (Monkhorst, Pack, 1976). Совокупность частот колебаний в Γ-точке не зависит от объема суперячейки, поскольку центр суперячейки остается неподвижным при заданных преобразованиях.

Вычисления проводились в два этапа: оптимизацию геометрических параметров решетки и затем, определение частот колебаний изотопологов. В процессе оптимизации геометрии по квази-ньютоновскому алгоритму определялись координаты ядер атомов и параметры решетки с минимальной энергией. Сходимость в процессе оптимизации оценивалась по среднеквадратичному отклонению (RMS) и абсолютному значению наибольшей компоненты градиентов и смещений ядер для всех атомов. Для максимальных и среднеквадратичных (RMS) градиентов выбраны значения (0.00012, 0.00003 a.u). Достоверность проведенных DFT вычислений подтверждается сравнением вычисленных в результате оптимизации геометрических параметров ячеек с экспериментальными.