ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КВАЗИСТАТИЧЕСКОГО НАГРУЖЕНИЯ ПОРИСТЫХ ФЛЮИДОНАСЫЩЕННЫХ ОБРАЗЦОВ ПОРОДЫ

Соловьев Сергей Александрович; Лисица Вадим Викторович

УДК 550.834

DOI: 10.33764/2618-981X-2021-2-2-298-311

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ КВАЗИСТАТИЧЕСКОГО НАГРУЖЕНИЯ ПОРИСТЫХ ФЛЮИДОНАСЫЩЕННЫХ ОБРАЗЦОВ ПОРОДЫ

Сергей Александрович Соловьев

Институт нефтегазовой геологии и геофизики им. А.А.Трофимука, Россия, 630090, г. Новосибирск, проспект Академика Коптюга, 3, кандидат физ.-мат. наук, старший научный сотрудник, e-mail: solovevsa@ipgg.sbras.ru

Вадим Викторович Лисица

Институт нефтегазовой геологии и геофизики им. А.А.Трофимука, Россия, 630090, г. Новосибирск, проспект Академика Коптюга, 3, доктор физ.-мат. наук, зав.лаб., e-mail: LisitsaVV@ipgg.sbras.ru

В данной статье представлен численный алгоритм моделирования низкочастотной нагрузки пороупругих материалов, заполненных жидкостью, и оценки эффективных частотно-зависимых соотношений деформаций и напряжений для таких сред. Алгоритм решает уравнение Био в квазистатическом состоянии в частотном пространстве. В результате для каждой временной частоты приходится решать систему линейных алгебраических уравнений. Мы используем прямой решатель, основанный на разложении $ LU $, чтобы разрешить SLAE. Предложенный алгоритм позволяет проводить оценку тензора жесткости в широком диапазоне частот.

Ключевые слова: уравнение Био, пороупругость, конечные разности

NUMERICAL WAVE FIELDS QUASISTATIC MODELING IN FLUID-FILLED POROELASTIC MEDIA

Sergey A. Solovyev

Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, 3, Akademika Koptyuga Ave., PhD, Senior Researcher, e-mail: solovevsa@ipgg.sbras.ru

Vadim V. Lisitsa

Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, 3, Akademika Koptyuga Ave., DSc, Head of laboratory, e-mail: LisitsaVV@ipgg.sbras.ru

This paper presents a numerical algorithm to simulate low-frequency loading of fluid-filled poroelastic materials and estimate the effective frequency-dependent strain-stress relations for such media. The algorithm solves Biot equation in quasi-static state in the frequency space. As a result a system of linear algebraic equations have to be solved for each temporal frequency. We use the direct solver, based on the $LU$ decomposition to resolve the SLAE. According to the presented numerical examples the suggested algorithm allows reconstructing the stiffness tensor within a wide Frequency range.

Keywords: Biot equation, poroelasticity, finite differences

Интенсивное развитие технологий захоронения С02 [13], [7], разработки геотермальной энергии [19], [11] ставит перед методами сейсмического мониторинга сложные задачи - оценку мобильности флюидов и транспортных характеристик

коллектора. Скорости сейсмических волн на низких частотах практически нечувствительны к изменениям в структуре коллектора, вызванных замещением флюидом, частичного химического растворения карбонатной матрицы породы и прочее. Однако частотно-зависимые эффекты могут быть зарегистрированы и потенциально интерпретированы. В частности, изменения в составе флюида или геометрии поро-вого пространства значительно влияют на потоки флюида, индуцированные волнами (WIFF). Эти потоки возникают из-за локальных градиентов давления при распространении сейсмических волн в трещиновато-пористых средах [15], [17]. Обычно рассматриваются WIFF из трещин во вмещающую среду (FB-WIFF) и перетоки между трещинами (FF-WIFF). FB-WIFF появляются, если распространяется низкочастотная волна. В этом случае период волны достаточно велик для образования потока даже в средах с довольно низкой проницаемостью. Интенсивность FB-WIFF определяется контрастом сжимаемости между вмещающей породой и материалом, заполняющим трещины [9], [17]. Распространение высокочастотных сигналов вызывает FF-WIFF, определяемую свойствами материала, заполняющего трещины, а также локальной связностью трещин [17], [9], [6]. К сожалению, теоретические исследования этого эффекта включают рассмотрение относительно простых моделей среды. Причем связность трещин учитывается только для пар разнонаправленных трещин [9]. Численное исследование явления также ограничено такими критериями связности трещин [17], [9], за исключением исследования [8], в котором авторы применяют статистическое моделирование сети трещин и оценивают полученную связность трещин. Одной из причин этого является отсутствие эффективного численного алгоритма для моделирования распространения сейсмических волн в пороупругих средах или для решения уравнения Био в квазистатическом состоянии.

Постановка задачи

Уравнение Био в квазистатической простановке.

Рассмотрим квазистатические уравнения Био, описывающие процессы диффузии в пороупругих средах, заполненных жидкостью, в низкочастотных режимах [2], [1]. Мы имеем дело с декартовыми координатами и ограничиваемся рассмотрением двухмерного случая, поэтому уравнения можно записать следующим образом:

д_ дх

_д_ дх

д дх

д дг

(К + 2/Л^ + + аМ V и ' дх и дг

- + -

к дх дг J

д + —

&

/

( ди диг ^ дг дх

/

(дих диг Л дг дх

д + —

дг

, ди^ _ лдиг лу(дых Л К— + (К + 2/)—- + аМ —- + —^ дх дг дх дг

аМ

(ди„ ди Л „ (^

дх дг

+М

дх дг

(ди„ ди„ Л „ (^

■ + -дх дг

+М

дх дг

■ Л

= т —

к

■ Л

= т—

к

= 0,

где й = (их,иг)т - вектор смещений точек матрицы, ^ = - вектор относи-

тельных смещений жидкости относительно матрица, Яи - коэффициент Ламе неосушенной породы, ¡- модуль сдвига, а - параметр Био-Уиллиса, 77 - динамическая вязкость жидкости, к- абсолютная проницаемость породы, с- временная частота. Значения параметров Ли, л, М и а обычно оцениваются по объемным модулям осушенной породы КЛ, флюида К1 и твердой матрицы К [12]:

К 2

а = 1--, Х = К — л, М = БК /а,

-ту ■'и. и 4 ' и '

з

Е_ 1/Кй -1/К к = К,

1/К -1/К +р(1/Кг -1/К) и 1 - Ба

Эффективная вязкоупругая среда

Для построения эффективных макроскопических частотно-зависимых упругих моделей необходимо потребовать, чтобы для любых напряжений, приложенных к единице объема, средние деформации в эффективной вязкоупругой и исходной пороупругих моделях совпадали. Квазистатическое приближение вязко-упругой среды имеет вид:

дх

—

дх

-у.

—у.

Си(с)-х + ^(с)-2

С55(с)

дх -г дх

-г

д + —

дг

С55(с)

( ду —у2 ^ -г дх

д + —

дг

—у -у

С,з(с) + Сзз (с )

дх дг

= 0,

(2)

где V = (ух,у2)г- вектор смещения, а тензор С - частотно-зависимый тензор жесткости. Уравнение определено в прямоугольной области Б = [1}х,ь2х]х[1}2,ь]]. В качестве граничных условий используются условия вида а-п = а0 на границе Ж). В

этих обозначениях И - внешняя нормаль, а а - тензор напряжений. Рассмотрим три базовых нагрузки: 1) Сжатие по направлению X:

д-„

—у.

О» = С„(с)—х + С,з(с)-^

дх -г

-у —у

1. агг = С,з(с)+ Сзз(с)

дх -г

Ъ = С»« + ^

I -г дх

—у„

—у.

= Ф., ^ = С„(с)—х + С,з(с)-^

дх -г

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

—у -у

= 0, агг = С,з(с)—+ Сзз(с)—^ дх -г

= Фх

= 0,

(3)

-в

1

х

2 х

2 г

2 = Ь

2=Ь

г

Решение задачи (2),(3) может быть построить аналитически, в случае постоянных коэффициентов: о^(х,г) = фх, о^(х,г) = 0, охг(х,г) = 0 . Таким образом, компоненты тензора деформаций можно представить в виде

Ехх = ЯцФх , Е1В = £13фх , = £15фх ,

(4)

где £ - компоненты тензора податливости (обратные к тензору жесткости).

Если начальные нагрузки известны и деформации вычислены, можно решить уравнения относительно компонентов тензора податливости. 2) Сжатие по направлению 7:

^ , ч дvv ^ . . дк

Охх = Си(т) + Си(®)

дх дг

д- дО = С13(т)^Т + Сзз(т)—^-

дх дг

х=Ь

г=Ь

= 0, о„ = с„(т) ^ + С1з(т)д

дх дг

д— д-= ф, О = С13(т)^Т + Сзз(т)—^

дх дг

= 0,

х=ь2

= Ф,

г=¿2

о. = С55(т)

\ дг дх

= 0,

дЭ

(5)

Решая задачу (2), (5) получаем:

е» = £13ф , = £33ф , = £35ф , (6)

Таким образом, можно восстановить второй столбец тензора податливости.

3) Касательные напряжения:

4)

^ , ч д-, _ . . д-у

Ох = Сп(т) + С^т)

дх дг

д- д-2. огг = С,3(т) д-*- + С33 (т) д.*.

дх дг

= 0, о„ = Сп(т)+ С,3(т)д дх дг

= 0, огг = С,3(т) % + С33 (т)д дх дг

= 0,

( д- д-Ох* = С55(0) +

дЭ

^ дг дх

Решая задачу (2), (7), получаем:

£хх = ^ = S35¥, =

(7)

3.

(8)

Таким образом, можно восстановить третий столбец тензора податливости.

Чтобы построить эффективную модель анизотропной вязкоупругой среды, необходимо решить три краевые задачи для системы (1) с граничными условиями (3), (5), и (7) соответственно. По решению каждой задачи нужно построить компоненты тензора деформаций, после чего решить системы уравнений (4), (6) и (8) относительно компонент тензора податливости. Однако, чтобы получить

х=ь

2=Ь

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

м>-п |ал=0

Для упрощения интерпретации результатов удобно рассматривать скорость и затухание сейсмических волн в эффективных вязкоупругих средах [4], [21], [23], а не компоненты тензора жесткости. Для этого необходимо разрешить дисперсионное соотношение для вязкоупругого волнового уравнения. Однако полученная модель является анизотропной, где скорость и затухание зависят от направления распространения. Поэтому мы ограничимся рассмотрением частных случаев, то есть скорость и затухание квазипродольной или qP-волны, распространяющейся в направлениях х и z, и квазипоперечной или qS-волны, распространяющейся в направлении х. В двумерных ортотропных средах скорости qS-волны в направлениях х и z совпадают. Для оценки скоростей и факторов качества воспользуемся формулами [3], [21]:

V = px Vpz = C C33 p ' V =

Qpx *Cn 3C,/ Qpz =, Qs KC55 3C55.

Наиболее трудоемкой частью предлагаемого подхода к осреднению является численное решение уравнения (1), которое описано ниже.

Конечно-разностная аппроксимация

Для аппроксимации уравнения (1) внутри области D = [L1,lLx]x[Lz,l2] мы предлагаем использовать равномерную прямоугольную сетку с шагами \ и Ъ2. Предположим, что границы области имеют полуцелые координаты; кроме того, L = Х/2, Ll = XNX -1/2, L = z1/2, Ll = zNt _1/2, где Nx и Nx - номера узлов сетки в соответствующем пространственном направлении. Схема расчетной области и сетки представлена на рис. 1.

Мы определяем сеточные функции на сдвинутых сетках по правилу:

(Ux \ j+1/2 = ux (xi , Zj+1/2 ) (wx )г, j+1/2 = wx (xi , Zj+1/2 ) (Uz )г +1/2, j = Uz (Xi+1/2 , Zj )

(w ),+u2, j = w (2, zj )• Таким образом, первое и третье уравнения из (1) аппроксимируются в точках (x, Zj+1/2), а второе и четвертое - в точках (Xг'+1/2, zj). Все коэффициенты хранятся в полуцелых точках, предполагая, что они постоянны в ячейке сетки. Чтобы сохранить второй порядок сходимости, необходимо вычислить модуль сдвига в целочисленных точках по правилу [22], [14], [10], [20]:

1 1

■ + -

¡, 1 41 ¡1+1/2,1+1/2 ¡/-1/2,1+1/2 ¡1+1/2,1-1/2 ¡-1/2,1-1/2 у

Рис.1 Эскиз расчетной области и сетки.

Ячейка сетки и положения компонентов решения показаны на рис. 2.

Рис.2 Ячейка сетки и положение компонент решения.

1

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

Вх и Ь =

■Л+1/2,/ .Л

I-1/2, /

к

и

-х

+ОК)

(10)

х, г

I

где / - любая достаточно гладкая функция, а индексы I и J могут быть целыми или полуцелыми.

Конечно-разностная аппроксимация уравнения (1) определяет систему линейных алгебраических уравнений (СЛАУ) Ах = Ь размером N = 2(N -1)N + 2(N -1)N, свойства которых являются предметом исследования.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Решение СЛАУ

Свойства построенной системы линейных уравнений следует исследовать для двух разных случаев. Первый - ненулевая частота, второй - с = 0; т.е. статическая нагрузка.

Если частота больше нуля; т.е. о 0, матрица комплексная, несимметричная. Более того, из-за использования граничных условий Неймана (3), или (5), или (7) матрица вырождена. Ядро дифференциального оператора (1), (3) состоит из трех векторов:

Г -х Л Гс ] Г ^ Л Г 01 Г ^ Л Г Сз 1

и 0 и -С3

0 0 0

У ^ ^ У 0 J У ^ ^ у 0 J У ^ ^ У 0 J

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

Если частота равна нулю, правые части третьего и четвертого уравнений в (1) становятся тривиальными. Таким образом, можно ввести новую переменную

V = Ч-ч> =

дх дг

и упростить систему (1)

д_ дх

д дх

д дг

(А + 2и)—~х + АИ — + аМУ V и ; дх и дг

д + —

дг

и

' дих диг л кдг дх ^

и

ди ди

У дг дх J

д + —

дг

А— + (А + 2и)~и^ + аМУ и дх У и ' дг

аМ

' дих +диг ^ У дх дг J

' дих диг л у дх дг J

+МУ

= 0,

где присутствуют только три независимые переменные. Таким образом, ранг матрицы системы линейных уравнений составляет 3 / 4Ы.

Решение системы линейных уравнений с сингулярной несимметричной комплексной матрицей является сложной задачей для итерационных методов [18]. На сходимость итерационных решателей сильно влияет выбор предобуслав-ливателя, который не является основной темой данного исследования. Кроме того, мы имеем дело с двумерными задачами, поэтому прямые методы могут

быть эффективно применены для решения СЛАУ. В частности, мы используем разреженный прямой решатель Intel MKL PARDISO, который очень эффективно оптимизирован для архитектур Intel.

Численные эксперименты

Сначала мы проверили алгоритм на однородной модели пороупругой среды со следующими параметрами:

М = 5.7 • 109 Pa, а = 0.87,п = 3 • 10-3Pa-s,K = 10-12т2Д0 = 6.09

_ in-12™2

kg

109 Pa, M = 6.72 • 109 Pa, p = 2650 ■kg, p = 1040 kg/m3, T = 1.5 и ф = 0.3.

m

Последние четыре параметра - это плотность породы, плотность жидкости, извилистость и пористость соответственно. Они не используются в квазистатической модели, но нужны для оценки скорости распространения и затухания волны в динамическом состоянии.

Для моделирования квазистатического нагружения пористого материала, заполненного жидкостью, мы рассмотрели квадратную область размером 1 м и с дискретизацией ^ = 0.002 т. Таким образом, размер задачи составил

= = 500 узлов. Мы выбрали пять частот 0, 0,1, 1, 10, 100, 1000 Гц, вычислили тензоры жесткости, а затем оценили скорость и затухание быстрой Р-волны, как показано на рисунке 3. Изначально модель была однородной и изотропной, поэтому скорость не зависела от направления распространения. Единственным физическим фактором, вызывающим затухание в однородных поро-упругих средах, является течение Био на высоких частотах. Однако его нельзя разрешить в квазистатическом состоянии. Таким образом, затухание численного решения равно нулю для всего диапазона частот. Скорость численного решения занижена, но разница составляет около 0,1 м/с; т.е. относительная погрешность составляет 3.5 • 10-5, что является приемлемым уровнем.

Рис.3 Сравнение численных оценок (красный) с аналитическим решением (синий). На левом рисунке представлена фазовая скорость, на правом рисунке -затухание

Трещиноватая среда

Чтобы проиллюстрировать применимость представленного подхода к оценке дисперсии и диссипации сейсмических волн, распространяющихся в неоднородных пороупругих средах, мы рассмотрели трещиновато-пористую среду. Мы использовали те же модели, что описаны в [16]. Мы рассмотрели два ортогональных набора трещин. Длина трещины в обоих наборах была зафиксирована на уровне 50 мм, а толщина - 2 мм. Мы создали несколько типов моделей трещин в зависимости от средней длины перколяции, используя метод имитационного отжига. Примеры моделей представлены на рисунке [4]. После этого мы заполнили модель свойствами материала. Вмещающая порода была слабопроницаемой, тогда как материал заполнения трещины был относительно гидравлически мягким, чтобы поддерживать поток жидкости. Описание модели представлено в таблице 1. Наличие неоднородности проницаемости в модели при распространении сейсмических возникают локальные перетоки флюида или перетоки, обусловленные распространением волн (от английского Wave-induced fluid flows) [15], [5], что в свою очередь приводит к диссипации сейсмической энергии.

Таблица 1

Свойства материалов.

Параметр Скелет Жидкость

Динамическая вязкость жидкости, Па 0.001 0.001

Проницаемость к0, м2 10А-15 5.5*10А-13

Плотность, кг / м3 2458 2458

Коэффициент Ламе 7.159*10А9 2.40*10А10

Модуль сдвига, Па 3.0969*10А10 1.14*10А10

Константа Био и Уиллиса 0.2962 0.6078

Коэффициент накопления жидкости М, Па 2.01*10А10 9.48*10А10

Мы использовали предложенный алгоритм для оценки частотно -зависимых тензоров жесткости для всех шести типов моделей в частотном диапазоне V 6 [0,..., 1000] Гц. Кроме того, мы непосредственно моделировали распространение волн в трещиновато-пористых средах, заполненных жидкостью, используя подход, описанный в [16]. Распространение волны моделировалось для импульса Риккера с центральной частотой 1000 Гц. На рисунках [5] и [6] мы приводим оценки скорости и затухания, полученные двумя методами, соответственно. Полученные оценки скоростей занижены, тогда как затухание хорошо согласуется. В целом полученные результаты показывают, что увеличение длины перколяции трещиноватой системы вызывает увеличение затухания из-за FB-WIFF.

Рис. 4. Модели трещинных сред с различной длиной перколяции (увеличивающейся сверху вниз).

Рис. 5. Скорости волны qP для моделей с разной перколяцией. Линии представляют квазистатические оценки, маркеры используются для оценок при моделировании распространения волн.

Рис. 6. Затухание волны qP для моделей с разной перколяцией.

Линии представляют квазистатические оценки, маркеры используются для оценок при моделировании распространения волн.

Оценки производительности

Во всех описанных примерах использовалась сетка с 500 узлами в обоих пространственных направлениях. Мы вычислили эффективные тензоры жесткости для набора из 22 частот от 0 до 1000 Гц. Разовый запуск алгоритма (одна модель, одна частота) включает следующие шаги:

1. Построение матрицы A и правых частей b;

2. PARDISO. Перестановка строк-столбцов матрицы. Это предварительный шаг для уменьшения объема памяти и времени, необходимого для факторизации матрицы системы.

3. PARDISO. Факторизация или LU разложение матрицы A.

4. PARDISO. Решение системы.

5. Постобработка: проверка относительной невязки (| b — Ax | / | b | < 10-12), а затем построение тензора жесткости C.

Следует отметить, что при изменении частоты необходимо корректировать только главную диагональ матрицы A. Таким образом, шаг (1) может быть выполнен для нулевой частоты, при следующих циклах алгоритма (при изменении частоты) меняются только диагональные элементы. Шаг переупорядочения (2) зависит только от позиции ненулевых элементов A и может быть выполнен только один раз, но применяется ко всем w. В результате весь алгоритм (шаги 1-5) применяется только к первой частоте в необработанном виде, тогда как сокращенная версия алгоритма (шаги 3-5) может применяться ко всем остальным частотам.

В наших расчетах использовался Intel (R) Xeon (R) CPU E5-2690 v2 @ 3,00 ГГц с 20 ядрами. Время расчета для одной модели и всех 22 частот составляет ~ 140s.

Профилирование времени 22 прогонов:

1. Аппроксимация задачи: ~ 3s.

2. PARDISO подготовительный шаг («реордеринг»): « 4s.

3. PARDISO шаг факторизации. 22 х ( « 3-4 с.)

4. PARDISO шаг решения. 22 х ( « 0.8 с.)

5. Постобработка: 22 х ( « 0.4 с.)

Для решения задачи размером 5002 потребовалось около 8 Гб оперативной памяти, то есть задачи такого размера можно решать на персональном компьютере или использовать GP-GPU для повышения производительности.

Заключение

Мы представили численный алгоритм для моделирования низкочастотного нагружения пороупругих материалов, заполненных жидкостью, и оценки эффективных частотно-зависимых соотношений деформаций и напряжений для таких сред. Алгоритм включает в себя решение уравнения Био в квазистатической постановке. Задача параболическая, поэтому ее удобно решать в частотной области. В результате система линейных алгебраических уравнений должна быть решена для каждой временной частоты. Мы используем прямой решатель, основанный на LU разложении. Предложенный алгоритм позволяет восстановить тензор жесткости в широком диапазоне частот [0, ..., 1000] Гц.

Исследование выполнено при поддержке Российского фонда фундаментальных исследований № 20-45-540004. Моделирование проводилось с использованием вычислительных ресурсов Сибирского суперкомпьютерного центра СО РАН.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Biot, M.A. Theory of propagation of elastic waves in a fluid-saturated porous solid. ii. higher frequency range // Journal of the Acoustical Society of America. - 1956. - 28. - P.179 - 191.

2. Biot, M.A.: Theory of propagation of elastic waves in fluid-saturated porous solid. i. low-frequency range // Journal of the Acoustical Society of America 28. - 1956. - P.168 - 178.

3. Carcione, J.M., Cavallini, F. A rheological model for anelastic anisotropic media with applications to seismic wave propagation // Geophys. J. Int. - 1994. - 119. - P.338 - 348.

4. Christensen R.M. Theory of viscoelasticity, an introduction. - Academic press New York and London. - 1971. - 249 p.

5. Germn Rubino J., Guarracino L., Mller T.M., Holliger K. Do seismic waves sense fracture connectivity? // Geophysical Research Letters. - 2013 - 40(4). - P.692- 696.

6. Guo J., Rubino J.G., Glubokovskikh S., Gurevich B. Effects of fracture intersections on seismic dispersion: theoretical predictions versus numerical simulations // Geophysical Prospecting. -2017. - 65(5). - P. 1264 - 1276.

7. Huang F., Bergmann P., Juhlin C., Ivandic M., Lth S., Ivanova A., Kempka T., Henninges J., Sopher D., Zhang F. The first post-injection seismic monitor survey at the ketzin pilot co2 storage site: results from time-lapse analysis // Geophysical Prospecting. - 2018. - 66(1). - P.62 - 84.

8. Hunziker J., Favino M., Caspari E., Quintal B., Rubino J.G., Krause R., Holliger K. Seismic attenuation and stiffness modulus dispersion in porous rocks containing stochastic fracture networks // Journal of Geophysical Research: Solid Earth. - 2018. - 123(1). - P.125 - 143.

9. Kong L., Gurevich B., Zhang Y., Wan Y. Effect of fracture fill on frequency-dependent an-isotropy of fractured porous rocks // Geophysical Prospecting. - 2017. - 65(6). - P.1649 - 1661.

10. Lisitsa V., Podgornova O., Tcheverda V. On the interface error analysis for finite difference wave simulation // Computational Geosciences. - 2010. - 14(4). - P.769 - 778.

11. Marty N.C.M., Hamm V., Castillo C., Thiry D., Kervvan C. Modelling water-rock interactions due to long-term cooled-brine reinjection in the dogger carbonate aquifer (paris basin) based on in-situ geothermal well data // Geothermics - 2020. - 88. - 101899.

12. Masson Y.J., Pride S.R., Nihei K.T. Finite difference modeling of biot's poroelastic equations at seismic frequencies. // Journal of Geophysical Research: Solid Earth. - 111(B10). - P.305.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

13. Menke H.P., Reynolds C.A., Andrew M.G., Pereira Nunes J.P., Bijeljic B., Blunt M.J.: 4d multi-scale imaging of reactive flow in carbonates: Assessing the impact of heterogeneity on dissolution regimes using streamlines at multiple length scales. // Chemical Geology. - 2018. - 481. - P.27

- 37.

14. Moczo P., Kristek J., Vavrycuk V., Archuleta R.J., Halada L. 3d heterogeneous staggered-grid finite-differece modeling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities. // Bulletin of the Seismological Society of America. - 2002. - 92(8). -P. 3042 - 3066.

15. Muller T.M., Gurevich B., Lebedev M. Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks - a review. // Geophysics. - 2010. - 75(5). - 75A147 -75A164

16. Novikov M.A., Lisitsa V.V., Bazaikin Y.V. Wave propagation in fractured-porous media with different percolation length of fracture systems // Lobachevskii Journal of Mathematics. - 2020.

- 41(8). - P.1533-1544

17. Rubino J.G., Muller T.M., Guarracino L., Milani M., Holliger K. Seismo-acoustic signatures of fracture connectivity // Journal of Geophysical Research: Solid Earth. - 2014. - 119(3). -P.2252 - 2271

18. Saad Y. Iterative Methods for Sparse Linear Systems. - SIAM. - 2003.

19. Salaun N., Toubiana H., Mitschler J.B., Gigou G., Carriere X., Maurer V., Richard A. Highresolution 3d seismic imaging and refined velocity model building improve the image of a deep geothermal reservoir in the upper rhine graben. // The Leading Edge. - 2020. - 39(12). - P.857 - 863

20. Samarskii A.A. The theory of difference schemes // Pure and Applied Mathematics. - v.240.

- CRC Press - 2001.

21. Vavrycuk V. Velocity, attenuation, and quality factor in anisotropic viscoelastic media: A perturbation approach // Geophysics. - 2008. - 73(5). - D63 - D73.

22. Vishnevsky D., Lisitsa V., Tcheverda V., Reshetova G. Numerical study of the interface errors of finite-difference simulations of seismic waves // Geophysics. - 2014. - 79(4). - P.T219 -T232

23. Zhu Y., Tsvankin I. Plane-wave propagation in attenuative transversely isotropic media // Geophysics. - 2006. - 71(2). - P. T17 - T30.

REFERENCES

1. Biot, M.A. Theory of propagation of elastic waves in a fluid-saturated porous solid. II. higher frequency range // Journal of the Acoustical Society of America. - 1956. - 28. - P.179 - 191.

2. Biot, M.A.: Theory of propagation of elastic waves in fluid-saturated porous solid. I. low-frequency range // Journal of the Acoustical Society of America. - 1956. - 28. - P.168 - 178.