Nature of anisotropic response of fluid saturated medium to surface seismic wave propagation

Введение

Горные породы при решении различных задач геодинамики, геофизики и геомеханики часто рассматриваются в приближении изотропного линейно-упругого материала с постоянными скоростями упругих волн. Однако хорошо известно, что породы, слагающие верхний слой земной коры, могут быть весьма неоднородными и содержать дефекты сплошности различного типа: разномасштабные трещины, каверны, границы блоков и т.д. Эти дефекты могут оказывать значительное влияние на физические и механические свойства пород, вызывая, в частности, анизотропию скоростей упругих волн [Виноградов и др., 1989; 1992; Егоркин, Егоркин, 1986; Chesnokov, Zatsepin, 1991; Nur, Simmons, 1969; Nur, 1971; Bonner, 1974; Lockner et al., 1977; Sayers, 2002; Stanchits et al., 2006; Browning et al., 2018]. Индуцированная дефектами анизотропия скоростей упругих волн в горных породах связана с наличием преимущественной ориентации дефектов, которая в свою очередь вызвана процессами формирования и деформирования породы в поле приложенных гравитационных и тектонических сил. Степень анизотропии зависит от ряда факторов, таких как плотность трещин, их ориентация, величина и вид напряженно-деформированного состояния породы. По мере изменения этих факторов может меняться отклик породы на прохождение сейсмических волн или упругих волн другой природы. Анизотропия скоростей упругих волн, индуцированная дефектами сплошности, также наблюдается вблизи крупных активных разломов и других сред с высоким уровнем трещиноватости [Crampin 1987; Leary et al., 1990; Miller, Savage 2001; Peng, Ben-Zion, 2004]. Учет индуцированной анизотропии поврежденных (трещиноватых) пород имеет важное значение для различных геофизических и геомеханических приложений, таких как сейсморазведка и геомеханическое моделирование резервуаров углеводородов.

В отличие от линейно-упругих материалов изменение объема при деформировании горных пород зависит от компонент девиатора напряжений. Так при воздействии негидростатической нагрузки трещиноватые породы могут демонстрировать значительную объемную деформацию еще до достижения предельной нагрузки [Lockner et al., 1992; Schmitt, Zoback 1992; Renard et al., 2019]. Согласно классической теории пороупругости Био [Biot, 1941], изменение порового давления флюида в трещиновато-пористой среде при недренированных условиях (когда отток флюида отсутствует или затруднен) определяется исключительно изменением среднего напряжения

$\Delta P_f=B\Delta {\sigma }_m$, (1)

где: $P_f$ - поровое давление; ${\sigma }_m$ - среднее напряжение в твердом скелете; $B=\beta M/K_u$ - коэффициент Скемптона, определяемый коэффициентом Био $\beta$, модулем Био $M$ и недренированным модулем объемного сжатия $K_u$. Многочисленные лабораторные экспериментальные исследования пороупругого поведения горных пород показали, что девиаторные напряжения могут оказывать существенное влияние на поровое давление [Wang, 1997; Lockner, Stanchits, 2002; Hamiel et al., 2005], а их связь определяется уравнением Скемптона [Skempton, 1954]:

 $\Delta P_f=B\left(\Delta {\sigma }_m+2\left(A-\frac{1}{3}\right)\Delta {\sigma }_d\right)$, (2)

где ${\sigma }_d$ - дифференциальное напряжение; А - дополнительный коэффициент Скемптона. Таким образом, величина дополнительного коэффициента Скемптона определяет чувствительность флюидонасыщенной породы к изменению девиаторных напряжений и возможность описания ее поведения в приближении линейного пороупругого тела.

Подземные воды, заключенные в коллекторах различного типа, выступают своеобразным индикатором напряженного состояния водовмещающих пород. В естественных условиях пороупругий отклик среды на экзогенные и эндогенные факторы приводит к соответствующим вариациям уровня водоносных горизонтов (порового давления). Это явление лежит в основе многих исследований, которые можно разделить на три направления. К первому направлению относится определение типа коллектора на основе гармонического анализа по гидрогеологическому отклику на приливы и изменение атмосферного давления [Rahi, Halihan, 2013]. Для безнапорного водоносного горизонта барометрический отклик является частотно-зависимым в отличие от этого же параметра для напорного горизонта. Согласно низко- и высокочастотным барометрическим откликам можно судить об ограниченности водоносного горизонта [Lai et al., 2013].

Второе направление посвящено оценке фильтрационных свойств флюидонасыщенных массивов и зон динамического влияния разломов по отклику на изменение атмосферного давления и приливные воздействия. Отклик порового давления в закрытой скважине на медленные изменения напряжений может рассматриваться в рамках квазистатической пороупругой теории [Барабанов и др., 1988; Копылова, Болдина, 2019; Hsieh et al., 1987; Burbey, 2010; Xue et al., 2016]. Анализ приливного отклика широко используется для оценки коллекторских свойств водонасыщенной среды [Cutillo, Bredehoeft, 2011; Doan et al., 2006] и направлен на исследование фазового сдвига между приливными волнами, выделенными в вариациях уровня водоносных горизонтов, и объемной деформации пласта.

Третье направление составляют исследования, направленные на изучение гидрогеологических эффектов, вызванных землетрясениями, и установленных не только в сейсмоактивных регионах, но и в платформенных условиях [Вартанян, 2019; Волейшо и др., 2007; Горбунова и др., 2015; Копылова, Болдина, 2019; 2020; 2021; Kitagawa et al., 2011; Wang, Manga, 2010]. Многолетние наблюдения за реакцией подземных вод на землетрясения позволили выделить основные типы гидрогеологических эффектов: предвестники, косейсмические и постсейсмические вариации уровня подземных вод [Киссин, 2015]. Предвестники прослеживаются в виде гидрогеологических, гидрогеохимических и температурных аномалий в наблюдательных скважинах и источниках перед землетрясениями [Копылова, Болдина, 2021; Вартанян, 2002; 2019; Kopylova, Boldina, 2021]. Постсейсмические эффекты могут быть связаны с формированием магистрального разрыва, сопровождаемого системой оперяющих трещин, необратимыми деформациями в дискретной блоково-иерархической геосреде, локальными изменениями структуры порово-трещинного пространства после прохождения сейсмических волн. Соответственно, косейсмические вариации уровня воды или порового давления в закрытой скважине вызваны прохождением сейсмических волн и могут быть как скачкообразными, так и постепенными.

В приближении линейной пороупругости предполагается, что уровень воды должен реагировать только на прохождение Р-волн и волн Релея, имеющих существенную объемную деформационную составляющую. Сдвиговые волны (S-волны) и волны Лява не вызывают значительных изменений уровня воды, поскольку их распространение сопровождается сдвиговой деформацией среды, напрямую не сжимающей или не расширяющей порово-трещинное пространство. Однако, как и в случае с лабораторными экспериментами, многочисленные скважинные измерения показали, что уровень воды (поровое давление) в измерительных скважинах реагирует на прохождение волн Лява [Wang et al., 2009]. При этом наблюдается устойчивая однозначная связь между вызванными волной объемными и дифференциальными деформациями и колебаниями порового давления, определяемая также уравнением Скемптона, записанном в деформациях [Shalev et al., 2016a; 2016b; Lutzky et al., 2020].

Наличие преимущественной ориентации в зонах динамического влияния активных разломов может приводить к азимутальной зависимости пороупругой реакции водонасыщенного коллектора на проходящие волны [Crampin, 1989; Winterstein, 1990]. В работе [Shalev et al., 2016b] введены два коэффициента $M$ и $N$, связывающие изменение порового давления с объемной и сдвиговой компонентами деформации:

 $d{P}_f=Md{\epsilon }_v+Nd{\epsilon }_d$, (3)

где ${\epsilon }_v$ - объемная деформация; ${\epsilon }_d$ - дифференциальная деформация. На основе длительных скважинных измерений установлено, что величина коэффициента $N$ зависит от степени трещиноватости породы коллектора, и может быть использована для ее оценки по сейсмологическим и гидрогеологическим данным. Новые результаты мониторинга порового давления в скважине, пройденной вблизи активной разломной зоны в районе резервуара сброса сточных вод Арбакл (Оклахома, США), представленные в работе [Barbour, Beeler, 2021], показали, что отношение коэффициентов $N/M$ уравнения Скемптона в деформациях (3), оцененное по вариациям порового давления в скважине при прохождении телесейсмических волн от удаленных землетрясений, имеет азимутальную зависимость. Этот факт указывает на анизотропию пороупругого отклика флюидонасыщенного коллектора на прохождение сейсмических волн, которая не описывается уравнением Скемптона (3).

Целью настоящей работы является развитие теоретических представлений о связи вариаций порового давления с вариациями объемных и сдвиговых деформаций, вызванных прохождением сейсмических волн, учитывающих как степень трещиноватости анизотропного коллектора, так и ее преимущественную ориентацию. Новая формулировка уравнения Скемптона базируется на ранее разработанной авторами анизотропной пороупругой модели деформирования трещиновато-пористой флюидонасыщенной среды [Lyakhovsky et al., 2022a], учитывающей преимущественную ориентировку трещиноватости пород в поле приложенных напряжений.

Теоретическая модель пороупругого отлика флюидонасыщенной трещиноватой среды

Для произвольного трехосного напряженно-деформированного состояния среды уравнение Скемптона (2) было переписано через октаэдрическое напряжение ${\tau }_{oct}$ в виде [Henkel, 1960; Henkel, Wade, 1966]

 $\Delta P_f=B\left(\Delta {\sigma }_m+\frac{3A-1}{\sqrt{2}}\Delta {\tau }_{oct}\right)$, (4)

где ${\tau }_{oct}=\frac{1}{3}\sqrt{{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_1-{\sigma }_3\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2}$ - октаэдрическое напряжение; ${\sigma }_i$ - главные напряжения. При $A=\frac{1}{3}$ уравнение (4) сводится к уравнению линейной пороупругости (1) с одним коэффициентом Скемптона. В случае $A<\frac{1}{3}$ вариации порового давления ниже вариаций, предсказываемых теорией Био, что неоднократно подтверждалось в лабораторных экспериментах [Scmitt, Zoback, 1992; Lockner, Byerlee, 1994; Lockner, Stanchits, 2002; Paterson, Wong, 2005].

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

 $\Delta P_f=B\Delta {\sigma }_m+A^{*}{}_{ij}\Delta {\tau }_{ij}$, (5)

где ${\tau }_{ij}={\sigma }_{ij}-{\sigma }_m{\delta }_{ij}$ - компоненты девиатора напряжений; $A^{*}{}_{ij}$ - компоненты дополнительного тензорного коэффициента Скемптона. Неизвестные скалярный и тензорный коэффициенты Скемптона могут быть найдены в общем виде из решения системы двенадцати уравнений:

 $\left\{\begin{array}{c}B={\left.\frac{\partial P_f}{\partial {\sigma }_m}\right|}_{{\tau }_{ij}=\text{const}}\\ A^{*}{}_{ij}={\left.\frac{\partial P_f}{\partial {\tau }_{ij}}\right|}_{{\sigma }_m=\text{const}}\end{array}\right.$. (6)

Сложность системы уравнений (6) заключается в том, что в общем случае невозможно получить ее полную математическую формулировку.

Как было сказано ранее, колебания порового давления, вызванные распространением сдвиговых волн и волн Лява, значительно усиливаются в сильно трещиноватых породах [Shalev et al., 2016b]. Для поверхностных волн такого типа приращение полной деформации в системе координат, связанной с направлением распространения волны ( $x_1$ - направление распространения волны, $x_2$ - направление ортогональное к направлению распространения волны), может быть представлено в виде суммы приращения объемной деформации $d{\epsilon }_v$ и приращения сдвиговой компоненты $d{\epsilon }_{12}$, связанной с тангенциальной компонентой смещения в волне. Тогда соотношение (5) в системе координат, связанной с направлением распространения волны, может быть записано в виде

 $d{P}_f=B_{\Phi }d{\sigma }_m+A_{\Phi }d {\tau }_{oct}$, (7)

где: приращения среднего и октаэдрического напряжений определяются через приращения соответствующих деформаций:

 $\begin{array}{c}d{\sigma }_m=\frac{\partial {\sigma }_m}{\partial {\epsilon }_v}d{\epsilon }_v+\frac{\partial {\sigma }_m}{\partial {\epsilon }_{12}}d{\epsilon }_{12}\\ d{\tau }_{oct}=\frac{\partial {\tau }_{oct}}{\partial {\epsilon }_v}d{\epsilon }_v+\frac{\partial {\tau }_{oct}}{\partial {\epsilon }_{12}}d{\epsilon }_{12}\end{array}.$ (8)

Коэффициенты $B_{\Phi }$ и $A_{\Phi }$ в этом случае определяются уравнениями

$\begin{array}{c}{B}_{\Phi }=\frac{\beta M}{\frac{\partial {\sigma }_m}{\partial {\epsilon }_v}-\frac{\partial {\sigma }_m}{\partial {\epsilon }_{12}}\frac{\frac{\partial {\tau }_{oct}}{\partial {\epsilon }_v}}{\frac{\partial {\tau }_{oct}}{\partial {\epsilon }_{12}}}},A_{\Phi }=\frac{\beta M}{\frac{\partial {\tau }_{oct}}{\partial {\epsilon }_v}-\frac{\partial {\tau }_{oct}}{\partial {\epsilon }_{12}}\frac{\frac{\partial {\sigma }_m}{\partial {\epsilon }_v}}{\frac{\partial {\sigma }_m}{\partial {\epsilon }_{12}}}}\end{array}$, (9)

где $\beta$, $M$ - коэффициент и модуль Био соответственно. Таким образом, для нахождения коэффициентов в модифицированном уравнении Скемптона (7) необходимо выбрать определяющие соотношения, связывающие компоненты тензора напряжений и деформаций и учитывающие тензорный характер развития поврежденности (трещиноватости) во флюидонасыщенной геосреде.

Ранее в работе авторов [Lyakhovsky et al., 2022a] была предложена модель нелинейного анизотропного пороупругого тела с тензорными поврежденностью и пористостью. Степень трещиноватости описывается с помощью тензора поврежденности второго ранга ${\Omega }_{ij}$, главные направления которого задают направления ортотропии материала, а главные значения определяют изменение эффективной площади поперечных сечений, перпендикулярных к каждой из осей ортотропии. Преимущество такого представления заключается в возможности определения компонент тензора для любого набора разноориентированных ансамблей микротрещин, полученного по томографическим или петрофизическим данным. Для описания эффекта направленного уплотнения вместо скалярного параметра пористости вводится тензорная величина ${\Psi }_{ij}$, след которой отражает пористость материала [Пантелеев и др., 2022; Lyakhovsky et al., 2022b]. Разработанная модель наследует основные черты ранее разработанных нелинейных моделей деформирования геосреды со скалярными параметрами поврежденности и пористости [Hamiel et al., 2004; 2005], а взаимодействие двух введенных тензорных параметров позволяет описать различные сценарии неупругого деформирования трещиноватой пористой среды [Lyakhovsky et al., 2022a].

Свободная энергия деформируемого трещиноватого (поврежденного) пористого материала, согласно этой модели, имеет вид:

$\begin{array}{c}F=\left[\frac{{\lambda }_0}{2}{I}_1^2+{\mu }_0{I}_2\right]+\left[\gamma {\xi }_0{I}_2^{\left(\Omega \right)}-\gamma I_1^{\left(\Omega \right)}\sqrt{I_2^{\left(\Omega \right)}}\right]+\\ +\frac{M}{2}{\left[\beta I_1-\zeta +{\Psi }_{ij}{\delta }_{ij}\right]}^2,\end{array}$ (10)

где: ${\lambda }_0$, ${\mu }_0$ - коэффициенты Ламе однородного (неповрежденного) материала; $\gamma$ - дополнительный упругий модуль, определяющий зависимость упругих свойств от вида напряженно-деформированного состояния и усиление нелинейности деформационного отклика с ростом поврежденности; $I_1,I_2$ - первый и второй инварианты тензора деформаций; $I_1^{\left(\Omega \right)},I_2^{\left(\Omega \right)}$ - первый и второй инварианты симметризованного тензора эффективных деформаций ${\epsilon }_{ij}^{\left(\Omega \right)}=\frac{1}{2}\left({\epsilon }_{ik}{\Omega }_{kj}+{\Omega }_{ik}{\epsilon }_{kj}\right)$; ${\xi }_0$ - материальный параметр, контролирующий переход от залечивания микротрещин к их росту; $\zeta$ - объемная доля флюида (отношение объема флюида к элементарному объему материала). Третий член в соотношении (10) с модулем Био $M$ и коэффициентом Био $\beta$ отличается от классической линейной пороупругости только членом ${\Psi }_{ij}{\delta }_{ij}$. В исходном недеформированном состоянии в условиях полного насыщения и равенства нулю порового давления флюид занимает все поровое пространство, т.е. $\zeta ={\Psi }_{ij}{\delta }_{ij}$. В общем случае изменение тензорной величины ${\Psi }_{ij}$ описывает необратимое изменение геометрии порового пространства, что при неизменности объемной доли флюида ( $\zeta =\text{const}$ ) приводит к увеличению или уменьшению порового давления.

Дифференцирование свободной энергии по компонентам тензора деформации приводит к выражению для компонент тензора напряжений:

 $\begin{array}{c}{\sigma }_{lm}=\left[{\lambda }_0{I}_1{\delta }_{lm}+2{\mu }_0{\epsilon }_{lm}\right]-\gamma {\Omega }_{lm}\sqrt{I_2^{\left(\Omega \right)}}+\\ +\gamma \left({\xi }_0-\frac{1}{2}\frac{I_1^{\left(\Omega \right)}}{\sqrt{I_2^{\left(\Omega \right)}}}\right)\frac{\partial I_2^{\left(\Omega \right)}}{\partial {\epsilon }_{lm}}+\\ +\beta M\left(\beta I_1-\zeta +{\Psi }_{ij}{\delta }_{ij}\right){\delta }_{lm},\end{array}$ (11)

где $\frac{\partial I_2^{\left(\Omega \right)}}{\partial {\epsilon }_{lm}}=\frac{1}{2}({\Omega }_{li}{\epsilon }_{ik}{\Omega }_{km}+{\Omega }_{mi}{\epsilon }_{ik}{\Omega }_{kl})+\frac{1}{2}({\Omega }_{lj}{\epsilon }_{mk}\times$ $\times {\Omega }_{kj}+{\Omega }_{mj}{\epsilon }_{lk}{\Omega }_{kj})$. Дифференцирование по объемной доли флюида $\zeta$ приводит к выражению для порового давления

 $P_f=M\left(-\beta I_1+\zeta -{\Psi }_{lm}{\delta }_{lm}\right)$. (12)

Из соотношения (12) видно, что на величину порового давления оказывает влияние упругая деформация твердого скелета ( $-\beta I_1$ ), объемная доля флюида $\zeta$ и необратимое изменение порового пространства, описываемое сверткой ${\Psi }_{lm}{\delta }_{lm}$.

Чтобы вычислить коэффициенты $B_{\Phi }$ и $A_{\Phi }$ в терминах упругих свойств, компонент тензора поврежденности и текущего состояния, необходимо вычислить все производные, входящие в (9) с использованием определяющего соотношения (11).

Матрица жесткости $K_{ijkl}$ представляет собой приращение компонент тензора напряжений ${\sigma }_{lm}$ при малом приращении компонент тензора деформаций ${\epsilon }_{ij}$:

 $\begin{array}{c}{K}_{ijnm}=\frac{\partial {\sigma }_{ij}}{\partial {\epsilon }_{nm}}=({\lambda }_0+{\beta }^{2}M){\delta }_{ij}{\delta }_{nm}+\\ +2{\mu }_0 \left({\delta }_{in}{\delta }_{jm}+{\delta }_{im}{\delta }_{jn}\right)-\frac{\gamma }{2}\frac{{\Omega }_{ij}}{\sqrt{I_2^{\Omega }}}\frac{\partial I_2^{\Omega }}{\partial {\epsilon }_{nm}}+\\ +\gamma \left({\xi }_0-\frac{1}{2}\frac{I_1^{\Omega }}{\sqrt{I_2^{\Omega }}}\right)\frac{{\partial }^2{I}_2^{\Omega }}{\partial {\epsilon }_{ij}\partial {\epsilon }_{nm}} -\\ -\frac{\gamma }{2}\left(\frac{\partial I_1^{\Omega }}{\partial {\epsilon }_{nm}}\frac{1}{\sqrt{I_2^{\Omega }}}-\frac{1}{2}\frac{I_1^{\Omega }}{I_2^{\Omega }\sqrt{I_2^{\Omega }}}\frac{\partial I_2^{\Omega }}{\partial {\epsilon }_{nm}}\right)\frac{\partial I_2^{\Omega }}{\partial {\epsilon }_{ij}}.\end{array}$ (13)

С использованием матрицы жесткости приращение среднего напряжения может быть найдено как

 $d{\sigma }_m=-\frac{1}{3}tr\left[\frac{\partial {\sigma }_{ij}}{\partial {\epsilon }_{nm}}d{\epsilon }_{nm}\right]=-\frac{1}{3}{K}_{ijnm}{\delta }_{ij}d{\epsilon }_{nm}$. (14)

Раскладывая приращение полной деформации в виде $d{\epsilon }_{nm}=-\frac{1}{3}d{\epsilon }_v{\delta }_{nm}+d{\epsilon }_{nm}^d$ получаем выражение для недренированного модуля объемного сжатия $K_u$, являющегося производной среднего напряжения по объемной деформации:

 $\frac{\partial {\sigma }_m}{\partial {\epsilon }_v}=K_u=\frac{1}{9}{K}_{ijnm}{\delta }_{ij}{\delta }_{nm}$. (15)

С учетом (13) недренированный модуль объемного сжатия может быть найден через инварианты $I_1^{\left(\Omega \right)},I_2^{\left(\Omega \right)}$ и их производные из соотношения

 $\begin{array}{c}{K}_u={\lambda }_0+{\beta }^{2}M+\frac{2}{3}{\mu }_0-\\ -\frac{2\gamma }{9}\left(\frac{{\Omega }_{nn}}{\sqrt{I_2^{\Omega }}}{\Omega }_{mi}{\epsilon }_{ik}{\Omega }_{km}-\left({\xi }_0-\frac{1}{2}\frac{I_1^{\Omega }}{\sqrt{I_2^{\Omega }}}\right)\times \right.\\ \times \left.{\Omega }_{ij}{\Omega }_{ij}-\frac{I_1^{\Omega }{\left({\Omega }_{mi}{\epsilon }_{ik}{\Omega }_{km}\right)}^2}{2I_2^{\Omega }\sqrt{I_2^{\Omega }}}\right),\end{array}$ (16)

которое при $\gamma =0$ редуцируется к хорошо известному соотношению линейной пороупругости $K_u={\lambda }_0+{\beta }^{2}M+\frac{2}{3}{\mu }_0$.

Производная среднего напряжения по сдвиговой компоненте деформации ${\epsilon }_{12}$ может быть найдена через матрицу жесткости как

 $\frac{\partial {\sigma }_m}{\partial {\epsilon }_{12}}=-\frac{1}{3}{K}_{ij12}{\delta }_{ij}$. (17)

Соответственно, полное выражение для этой производной имеет вид:

 $\begin{array}{c}\frac{\partial {\sigma }_m}{\partial {\epsilon }_{12}}=-\frac{\gamma }{3}\left[-\frac{1}{2}\frac{{\Omega }_{kk}}{\sqrt{I_2^{\Omega }}}\frac{\partial I_2^{\Omega }}{\partial {\epsilon }_{12}}+4\left({\xi }_0-\frac{1}{2}\frac{I_1^{\Omega }}{\sqrt{I_2^{\Omega }}}\right)\times \right.\\ \times {\Omega }_{1m}{\Omega }_{m2}-\left({\Omega }_{mi}{\epsilon }_{ik}{\Omega }_{km}\right)\times \\ \times \left.\left(\frac{2{\Omega }_{12}}{\sqrt{I_2^{\Omega }}}-\frac{1}{2}\frac{I_1^{\Omega }}{I_2^{\Omega }\sqrt{I_2^{\Omega }}}\frac{\partial I_2^{\Omega }}{\partial {\epsilon }_{12}}\right)\right],\end{array}$ (18)

где $\frac{\partial I_2^{\Omega }}{\partial {\epsilon }_{12}}=\frac{1}{2}({\Omega }_{1i}{\epsilon }_{ik}{\Omega }_{k2}+{\Omega }_{i2}{\epsilon }_{ik}{\Omega }_{1k})+\frac{1}{2}({\Omega }_{1j}{\epsilon }_{k2}{\Omega }_{kj}+$ $+{\Omega }_{j2}{\epsilon }_{1k}{\Omega }_{kj})$.

Для вычисления производных октаэдрического напряжения перепишем его в произвольной системе координат, необязательно связанной с главными осями тензора напряжений:

$\begin{array}{c}{\tau }_{oct}=\frac{1}{3}\sqrt{{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_1-{\sigma }_3\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2}=\\ =\frac{1}{3}\sqrt{3S_2-S_1^2},\end{array}$ (19)

где $S_1={\sigma }_{kk}=-3{\sigma }_m$, $S_2={\sigma }_{ij}{\sigma }_{ij}$. Используя соотношения, приведенные выше, получаем:

 $\frac{\partial {\tau }_{oct}}{\partial {\epsilon }_v}=\frac{1}{18{\tau }_{oct}}\left[3\frac{\partial S_2}{\partial {\epsilon }_v}-2S_1\frac{\partial S_1}{\partial {\epsilon }_v}\right]=\frac{-{\tau }_{ij}{K}_{ijnm}{\delta }_{nm}}{9{\tau }_{oct}}$, (20)

 $\frac{\partial {\tau }_{oct}}{\partial {\epsilon }_{12}}=\frac{1}{18{\tau }_{oct}}\left[3\frac{\partial S_2}{\partial {\epsilon }_{12}}-2S_1\frac{\partial S_1}{\partial {\epsilon }_{12}}\right]=\frac{{\tau }_{ij}{K}_{ij12}}{3{\tau }_{oct}}$, (21)

где ${\tau }_{ij}$ - компонент девиатора напряжений.

Подстановка найденных производных среднего и октаэдрического напряжений по объемной деформации и сдвиговой компоненте приводит к выражениям для коэффициентов $B_{\Phi }$ и $A_{\Phi }$ в модифицированном уравнении Скемптона (7):

 $\begin{array}{c}{B}_{\Phi }=\frac{9\beta M{\tau }_{ij}{K}_{ij12} }{9K_u{\tau }_{ij}{K}_{ij12}-K_{ij12}{\delta }_{ij}{\tau }_{ij}{K}_{ijnm}{\delta }_{nm}}\\ A_{\Phi }=\frac{9\beta M{\tau }_{oct}{K}_{ij12}{\delta }_{ij}}{9K_u{\tau }_{ij}{K}_{ij12}-K_{ij12}{\delta }_{ij}{\tau }_{ij}{K}_{ijnm}{\delta }_{nm}}\end{array}$. (22)

Видно, что найденные коэффициенты зависят от текущего уровня приложенных напряжений ${\tau }_{ij}$ и ${\tau }_{oct}$, а не только от физических (материальных) параметров, что крайне неудобно при использовании их для анализа конкретных геофизических данных. От этого недостатка можно избавиться, перейдя от уравнения Скемптона в напряжениях (7) к уравнению Скемптона в деформациях

 $d{P}_f=M_{\Phi }d{\epsilon }_v+N_{\Phi }d {\epsilon }_{12}$. (23)

Коэффициенты $M_{\Phi }$ и $N_{\Phi }$ могут быть найдены через коэффициенты уравнения (7) с использованием соотношений

 $M_{\Phi }=B_{\Phi }{K}_u, N_{\Phi }=A_{\Phi }\frac{\partial {\tau }_{oct}}{\partial {\epsilon }_{12}}=A_{\Phi }\frac{{\tau }_{ij}{K}_{ij12}}{3{\tau }_{oct}}.$ (24)

Следуя работе [Barbour, Beeler, 2021], найдем отношение этих коэффициентов:

 $\begin{array}{c}\frac{N_{\Phi }}{B_{\Phi }{K}_u}=\frac{1}{3K_u}{K}_{ij12}{\delta }_{ij}=\\ =\frac{1}{3K_u}\left(K_{1112}+K_{2212}+K_{3312}\right).\end{array}$ (25)

Полученное выражение для $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ обладает рядом преимуществ. Во-первых, оно определяется только материальными параметрами, величиной главных компонент и ориентацией тензора поврежденности, а также видом напряженно-деформированного состояния. Во-вторых, полученное выражение имеет прозрачный физический смысл. Отношение коэффициентов, определяющих пороупругий отклик флюидонасыщенного поврежденного материала на изменение объемной и сдвиговой деформации, равно отношению среднего модуля, связывающего диагональные компоненты тензора напряжений со сдвиговой компонентой деформации ${\epsilon }_{12}$, к недренированному модулю объемного сжатия. В случае $\frac{N_{\Phi }}{B_{\Phi }{K}_u}\ll 1$ можно говорить о нечувствительности поврежденной флюидонасыщенной среды к прохождению сдвиговых волн и волн Лява, а при $\frac{N_{\Phi }}{B_{\Phi }{K}_u}\approx 1$ о высокой чувствительности среды, сопоставимой с пороупругим откликом на прохождение Р-волн и волн Релея.

Численное моделирование пороупругого отклика водонасыщенной трещиноватой среды в окрестности активного разлома

Для оценки влияния трещиноватости, имеющей преимущественную ориентацию, на отношение коэффициентов $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ в модифицированном уравнении Скемптона (23) в случае воздействия на трещиноватую флюидонасыщенную среду в районе активного разлома поверхностной волны рассмотрим следующую задачу. Сдвиговая разломная зона находится в поле природных напряжений, ориентированных следующим образом (рис. 1): вертикальная компонента тензора напряжения (отрицательная ${\sigma }_v$ ) действует вдоль оси Z, минимальное горизонтальное напряжение $S{h}_{\text{min}}={\sigma }_v+\Delta$ ориентировано вдоль оси X, максимальное горизонтальное напряжение $S{h}_{\text{max}}={\sigma }_v-\Delta$ - вдоль оси Y. Параметр $\Delta$ определяет уровень напряжений, обеспечивающих выполнение критерия активизации разломной зоны. Разломная зона ориентирована к направлению действия максимального горизонтального напряжения под углом Кулона-Мора $\theta$, определяемым углом внутреннего трения среды $\phi$ как $\theta =\frac{\pi }{4}\pm \frac{\phi }{2}$. Направление прихода поверхностной волны задается углом $\omega \in \left[\mathrm{0,2}\pi \right]$ и определяет локальную систему координат $O{X}^{*}{Y}^{*}$, в которой будет оцениваться отношение коэффициентов $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$, что влечет необходимость вращения всех тензорных величин на соответствующий угол.

Рис. 1. Ориентация зоны сдвига и два варианта ориентации трещиноватости: А – изотропная трещиноватость (равновероятно ориентированная); В – трещиноватость, ориентированная параллельно оси максимального горизонтального сжатия.

Для нахождения зависимости $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ от угла прихода волны $\omega$ были выбраны следующие значения материальных параметров, характерные для песчаников [Wang, 2000]: параметр Ламе ${\lambda }_0$ = 30 МПа, отношение параметров Ламе μ0 = $=2.5{\lambda }_0$, модуль Био $M=2.84{\lambda }_0$ и коэффициент Био $\beta =0.7$, угол внутреннего трения $\phi =40°$. Для принятых значений материальных параметров получаем ${\xi }_0=-1$ и дополнительный модуль близок модулю сдвига $\gamma ~{\mu }_0$.

Будем рассматривать два случая ориентации трещиноватости (поврежденности) (см. рис. 1): изотропную трещиноватость, характеризующуюся равновероятным распределением микротрещин по ориентации (случай А), трещиноватость, ориентированную параллельно оси максимального горизонтального сжатия (случай В). Необходимо отметить, что второй случай ориентации трещиноватости соответствует многочисленным экспериментальным наблюдениям [Reches, Lockner, 1994] и ранее был теоретически предсказан в рамках нелинейной модели упругости с тензорным параметром поврежденности как оптимальный угол ориентации трещиноватости в районе сдвигового разлома, доставляющий максимум скорости диссипации энергии [Пантелеев, Ляховский, 2022]. Оба случая характеризуются заданным средним уровнем поврежденности $\alpha =\frac{1}{3}\left({\Omega }_x^2+{\Omega }_y^2+{\Omega }_z^2\right)=0.3$. Соответственно, компоненты тензора поврежденности для случая изотропной поврежденности А равны ${\Omega }_x={\Omega }_y={\Omega }_z=\sqrt{0.3}$, для случая B - ${\Omega }_y={\Omega }_z=\mathrm{0,} {\Omega }_x=\sqrt{0.9}$. Реальные зоны разломов, вероятно, будут иметь сложные распределения ориентации трещиноватости, однако и в этом случае можно выделить преимущественную ориентацию. Несмотря на одинаковую среднюю поврежденность для случаев А и В, величина параметра $\Delta$, определяющая условие начала роста поврежденности, будет различной: $\frac{\Delta }{{\sigma }_v}=0.47$ для случая А, $\frac{\Delta }{{\sigma }_v}=0.39$ для случая В.

Процедура расчета зависимости отношения коэффициентов $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ от угла прихода волны $\omega$ включает несколько этапов. Сначала для заданных главных напряжений, ориентации разломной зоны, ориентации трещиноватости и компонент тензора поврежденности в предположении тонкого поврежденного слоя рассчитываются компоненты тензора деформации в системе координат разломной зоны. Подробное описание алгоритма нахождения компонент тензора деформаций в тонком поврежденном слое представлено в работе [Пантелеев, Ляховский, 2022]. Далее для заданного угла прихода $\omega$ рассчитываются все компоненты тензора поврежденности и деформации в системе координат, связанной с направлением прихода волны (углом $\omega$ ). После этого, согласно соотношениям (13) и (16), вычисляются компоненты матрицы жесткости $K_{ijnm}$ и недренированный модуль объемного сжатия $K_u$, которые подставляются в соотношение (25) для расчета $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$. Процедура повторяется для следующего значения угла прихода волны $\omega$.

На рис. 2 приведена зависимость отношения коэффициентов модифицированного уравнения Скемптона $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ от угла прихода волны $\omega$. В отличие от того, что предсказывает традиционное уравнение Скемптона (соотношение (7)), представленная зависимость имеет периодический характер, симметричный относительно нуля. Относительные амплитуды отклика порового давления по модулю при варьировании угла прихода волны представлены на рисунке 2б штрихпунктирной линией. Так, существуют направления (на рис. 2б обозначены $n_k$ ), для которых отклик поврежденной среды на прохождение поверхностной волны будет полностью отсутствовать. Эти направления соответствуют направлениям действия главных горизонтальных напряжений. Наряду с этим есть направления $m_k$, для которых такой отклик будет максимальным.

Рис. 2. Зависимость отношения $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ от угла прихода волны для случая изотропной поврежденности (а), углы максимального по модулю ($m_i)$ и нулевого ( $n_i$ ) отклика порового давления на прохождение поверхностной волны (б).

Эти направления соответствуют углам $\frac{\pi }{4}\pm \frac{\pi }{2}$, т.е. направлениям действия максимальных касательных напряжений в среде. При этом знак отклика, определяющий увеличение или уменьшение порового давления при прохождении сейсмической волны, зависит от того, совпадает ли направление сдвига в волне с направлением сдвига в разломной зоне.

Также необходимо отметить, что коэффициент $N_{\Phi }$ более чем на порядок величины меньше коэффициента $B_{\Phi }{K}_u$, определяющего реакцию порового давления на изменение объемной деформации.

В случае анизотропной поврежденности (случая В) зависимость отношения коэффициентов $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ от угла прихода волны $\omega$ имеет такую же форму (рис. 3) и углы нулевого и максимального откликов, как и в случае А. При этом максимальная амплитуда отклика порового давления фактически на порядок превышает амплитуду для случая изотропной поврежденности (случай А).

Рис. 3. Зависимость отношения $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ от угла прихода волны для случая трещиноватости, ориентированной параллельно оси максимального горизонтального напряжения, случай B (а), углы максимального по модулю ($m_i)$ и нулевого ( $n_i$ ) отклика порового давления на прохождение поверхностной волны (б).

Увеличение амплитуды отклика для случая В связано с тем, что данное направление трещиноватости (параллельное направлению действия максимального горизонтального напряжения) является оптимальным с точки зрения роста трещин, согласно работе [Reches, Lockner, 1994], которые в дальнейшем объединяются в более крупные трещины вдоль простирания разлома, накапливая сдвиг. Необходимо отметить, что совпадение углов нулевого и максимального по модулю откликов со случаем изотропной поврежденности обусловлено в данном случае соосностью тензора действующих напряжений и тензора поврежденности. Для преимущественной трещиноватости, повернутой относительной направления действия главных напряжений, зависимость отношения коэффициентов $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ от угла прихода волны $\omega$ будет несимметрична относительно нуля и будет иметь сдвиг по фазе (нулевой отклик будет наблюдаться для углов, отличных от $0\pm \frac{\pi }{2}n$ ).

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

Верификация модифицированного уравнения Скемптона на данных скважинного мониторинга

Для верификации модифицированного уравнения Скемптона, у которого коэффициенты зависят как от вида напряженно-деформированного состояния, так и от величины главных компонент и ориентации тензора поврежденности, были использованы данные комплексного скважинного мониторинга резервуара для утилизации сточных вод Арбакл, Оклахома, США [Barbour, Beeler, 2021]. На основе анализа сейсмологических, тектонофизических и геологических данных в работе [Kolawole et al., 2019] показано, что разломы в Оклахоме простираются с северо-востока на юго-запад и с северо-запада на юго-восток. На рис. 4 приведена сводная схема возможных ориентаций разломной зоны, ориентаций направления действия максимального горизонтального напряжения и расположения эпицентров сейсмических событий относительно наблюдательной скважины в районе резервуара Арбакл. Видно, что зоны возможного расположения разломов благоприятно ориентированы относительно направления максимального горизонтального напряжения (под острым углом $\theta$ ) [Alt, Zoback, 2017]. Сейсмологические данные для указанных на рис. 4 сейсмических событий были использованы для оценки отношения коэффициентов $\frac{M}{N}$ уравнения Скемптона в деформациях (соотношение (3)). В результате проведенного анализа установлено, что величина $\frac{M}{N}$ имеет явно выраженную азимутальную зависимость [Barbour, Beeler, 2021]. Верификация предложенного модифицированного уравнения Скемптона заключалась в сопоставлении теоретической зависимости отношения коэффициентов от угла прихода волны, построенной с использованием соотношения (25), с величинами, оцененными по данным скважинного мониторинга вариации порового давления при прохождении сейсмических волн от семи выделенных событий и представленными в таблице № 2 в работе [Barbour, Beeler, 2021].

Рис. 4. Интервалы возможной ориентации разломов, ориентации оси максимального горизонтального напряжения и положение сейсмических событий относительно наблюдательной скважины в районе резервуара Арбакл [Barbour, Beeler, 2021].

Для построения теоретической зависимости отношения коэффициентов $B_{\Phi }{K}_u$ от угла прихода волны по соотношению (25) необходимо иметь информацию о трех углах и двух скалярных величинах: угле $\omega$ между направлением $S{H}_{\text{max}}$ и направлением прихода сейсмической волны, угле ориентации разлома $\theta$ относительно направления $S{H}_{\text{max}}$, ориентации трещиноватости $\varphi$ относительно направления $S{H}_{\text{max}}$, среднем уровне поврежденности $\alpha$ и отношении $\frac{\Delta }{{\sigma }_v}$. В связи с тем, что максимальная величина $\frac{M}{N}$, оцененная для выделенных сейсмических событий, близка к единице, то в качестве первого приближения был выбран случай В с трещиноватостью, ориентированной параллельно направлению $S{H}_{\text{max}}$. При этом было учтено, что по оценкам, приведенным в работе [Alt, Zoback, 2017], в районе резервуара Арбакл имеет место следующее соотношение природных напряжений ${\sigma }_v=S{H}_{\max }>S{H}_{\min }$, ${\sigma }_v=10$ МПа. Используя это соотношение, получаем безразмерное соотношение $\frac{S{H}_{\min }}{{\sigma }_v}=1-\frac{\Delta }{{\sigma }_v}$, зависящее только от параметра $\Delta$. Оценка неизвестных параметров осуществлялась путем решения задачи минимизации невязки между рассчитанной по соотношению (25) теоретической кривой и оценками отношения $\frac{M}{N}$ для семи представленных сейсмических событий:

 $\begin{array}{c}{\displaystyle \sum }_{i=1}^7{\left(\frac{N_{\Phi }}{B_{\Phi }{K}_u}\left(\omega ,\Delta ,\theta ,\alpha ,\varphi \right)-{\left(\frac{M}{N}\right)}_i\right)}^2\to \min ,\\ \begin{array}{c}0\le \omega \le 2\pi , 0\le \Delta \le 5 \text{ÌÏà}, \\ 24°\le \theta \le 37.5°, \\ 0.01\le \alpha \le \mathrm{0.8,}- 10°\le \varphi \le 10°.\end{array}\end{array}$ (26)

При решении задачи минимизации методом прямого поиска (методом Хука-Дживса) учитывались следующие ограничения на величину управляющих параметров:

• угол $\omega$ прихода волны варьируется в интервале от 0 до $2\pi$;
• параметр $\Delta$ может принимать значения от 0 до 5 МПа;
• угол наклона разломной зоны относительно направления $S{H}_{\text{max}}$ $\theta$ изменяется в интервале от 24 до 37.5°, что, с одной стороны, определяется имеющимися оценками [Kolawole et al., 2019] (см. рис. 4), а с другой стороны определяет диапазон вариации угла внутреннего трения среды от 15 до 42°;
• средняя поврежденность $\alpha$ может принимать значения в интервале от 0.01 до 0.8;
• угол наклона трещиноватости $\varphi$ относительно направления $S{H}_{\text{max}}$ варьируется в пределах ± $10°$.

В результате решения задачи минимизации получено, что наилучшее соответствие между теоретической кривой и оцененными по данным скважинного мониторинга значениями отношения коэффициентов уравнения Скемптона имеет место для следующих величин заданных (управляющих) параметров: параметр $\Delta$ = 3.88 МПа, ориентация разломной зоны $\theta =$ 32° относительно направления $S{H}_{\text{max}}$, что соответствует углу внутреннего трения в 26°, средний уровень поврежденности $\alpha$ = 0.4, трещиноватость ориентирована относительно направления $S{H}_{\text{max}}$ под углом $\varphi$ = 6°. Таким образом, полученный набор ориентаций показывает, что ситуация в районе резервуара Арбакл действительна близка к случаю В, рассмотренному ранее.

На рис. 5 представлено сравнение найденной теоретической зависимости $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ от угла прихода сейсмической волны $\omega$, рассчитанной для условия резервуара Арбакл, и оценок, полученных по сейсмологическим данным для семи сейсмических событий [Barbour, Beeler, 2021]. Видно, что отношение коэффициентов в уравнении Скемптона для шести из семи сейсмических событий с погрешностью менее 3% описывается построенной теоретической кривой. При этом для седьмого события имеет место ошибка в исходной оценке $\frac{M}{N}$, так как событие с моментной магнитудой 7.7 имеет близкий азимут, но существенно большее отношение коэффициентов $\frac{M}{N}$. Таким образом, полученный результат показывает, что модифицированное выражение для коэффициентов в уравнении Скемптона, построенное на основе нелинейной модели анизотропной пороупругости с тензорной поврежденностью, позволяет описать азимутальную зависимость отклика порового давления на прохождение сейсмических волн, обнаруженную в результате анализа данных мониторинга порового давления в реальном объекте. Предложенное соотношение (25) может быть использовано для оценки ориентации трещиноватости флюидонасыщенного коллектора в районе активных разломных зон по данным вариаций порового давления на прохождение поверхностных сейсмических волн.

Рис. 5. Сравнение расчетной кривой (соотношение (25)) и сейсмологических данных для отклика порового давления ( $\frac{N_{\Phi }}{B_{\Phi }{K}_u}$ ), индуцированного преимущественно ориентированной трещиноватостью пород.

Заключение

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

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

Финансирование работы

Работа выполнена при финансовой поддержке Российского научного фонда в рамках научного проекта № 19-77-30008.