Моделирование направления осей главных напряжений горных пород при подготовке землетрясений

Полный текст

Введение

Процесс изменения напряженно-деформированного состояния горных приводит к возникновению акустического излучения, которое называется акустической эмиссией горных пород или геоакустической эмиссией. Источники этого излучения имеют дислокационную природу и могут быть обусловлены упругими [1], пластическими [2] или гранулированными свойствами среды [3].

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

Одним из видов аномалий, проявляющихся на заключительной стадии подготовки землетрясений, является возникновение ярко выраженных направлений акустической активности [7]. Основная гипотеза возникновения этого явления заключается в том, что влияние очага готовящегося землетрясения приводит к формированию постоянного направления осей главных напряжений в точке наблюдений. Направления этих осей в свою очередь определяют преимущественную ориентацию источников акустического излучения.

В настоящей работе для подтверждения выдвинутой гипотезы произведено математическое моделирование направления осей главных напряжений, возникающих в горных породах при подготовке землетрясений вблизи полуострова Камчатка. Результаты моделирования сопоставлены с экспериментальными оценками направления осей главных напряжений по данным аномалий направленности акустической эмиссии горных пород в пункте наблюдений «Микижа» (52,99° с. ш., 158,22° в. д.) [8].

Физическая постановка задачи

Направленность геоакустической эмиссии. В процессе генерации акустического излучения горными породами преобладают сдвиговые источники. Это связано с тем, что прочность пород по отношению к касательным напряжениям меньше, чем к сжатию. Преимущественная ориентация таких источников определяется направлением максимальных касательных напряжений, ориентированных в основном под углом 45° к оси наибольшего сжатия [9, 10]. При этом максимумы продольных акустических колебаний соответствуют направлениям, перпендикулярным этой оси для волны сжатия и совпадающим с ней для волны разряжения. Диаграммы направленности акустических колебаний сдвиговых источников схематично представлены на рис. 1.

 

Рис. 1. Диаграммы направленности продольных (a) и поперечных (b) акустических колебаний сдвиговых источников. Стрелками показана ось наибольшего сжатия.

[Figure 1. Directivity diagrams of longitudinal (a) and transverse (b) acoustic oscillations of shear sources. Arrows indicate the strongest compression axis.]

 

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

Так, например, проведенный анализ аномалий направленности геоакустической эмиссии, зарегистрированной в пункте наблюдения «Микижа» (52,99° с. ш., 158,22° в. д.) на полуострове Камчатка в2008 $–$ 2012 гг., показал следующие результаты [8].

Максимумы акустической активности регистрировались с направлений близких к 40 и 220, несколько реже $–$ с направлений 140 и более широкого диапазона 300 $–$ 340. На основании этих данных были получены оценки направления осей главных напряжений. Диапазон направления осей главных напряжений изменялся в интервале углов 110 $–$ 140 со средним значением в 125 (рис. 2a). Полученный усредненный результат оценки ориентации оси совпал с генеральным направлением действия тектонических напряжений у побережья Южной Камчатки (рис. 2b).

 

Рис. 2. Оценка ориентации оси максимального сжатия перед 74 землетрясениями в период 2008 $–$ 2012 гг., ■$—$ область ориентации оси сжатия (a). Ориентация тектонических напряжений у побережья Камчатки, 🡪 🡨 $—$ направление осей горизонтального сжатия (b) [8].

[Figure 2. Estimate of the compression axis orientation before 74 earthquakes within the period 2008$—$2012 (a), ■$—$ the compression axis orientation region. Orientation of tectonic stresses near the Kamchatka peninsula, 🡪 🡨$—$ direction of horizontal compression axes (b) [8].]

 

Модель очага землетрясения. Очаг тектонического землетрясения представляет собой разрыв сплошности материала Земли, который возникает под действием упругих сдвиговых напряжений, накопленных в процессе тектонической деформации. В момент землетрясения происходит полное или частичное снятие накопленных напряжений в его очаге. Соответствующий очагу тектонического землетрясения разрыв является разрывом скольжения, то есть взаимное перемещение его берегов по нормали к поверхности разрыва равно нулю. Данное определение тектонического землетрясения было сформулировано Б. В. Костровым и основано на теории упругой отдачи Г. Ф. Рейда, которая гласит, что тектоническое землетрясение состоит в превращении части высвободившейся упругой потенциальной энергии в кинетическую энергию сейсмических волн.

Таким образом, деформации, возникающие при подготовке землетрясения, обусловлены приращением потенциальной энергией упругих деформаций ΔW, вызванным процессом подготовки землетрясения. Эта энергия больше, чем высвободившаяся сейсмическая энергия E. Величина η, равная отношению этих энергий, определяет эффективность снятия потенциальной энергии упругих деформаций и называется КПД землетрясения

$\eta \mathrm{<mo>=</mo>}\frac{E}{\Delta W}\mathrm{.}$ (1)

В механике сплошных сред процесс потери устойчивости среды изучается при неизменных свойствах материала среды и без нарушения его сплошности. В рамках этих ограничений описать очаг землетрясения можно через некоторую систему сил, распределенную по поверхности разрыва. Для описания произвольно ориентированного разрыва смещений в изотропной среде используют систему, состоящую из девяти пар двойных сил [11]. Такая система, схематически изображенная на рис. 3, соответствует модели очага землетрясения в виде двойного диполя. Пара двойных сил с моментами является минимально-необходимой системой сил для построения эквивалентного источника.

 

Рис. 3. Девять пар сил, необходимых для получения силового эквивалента при произвольно ориентированном разрыве смещений в среде.

[Figure 3. Nine pairs of forces required to obtain a force equivalent for an arbitrarily oriented displacement gap in the medium.]

 

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

Математическая постановка задачи

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

$\mu u_{i\mathrm{,}jj}+\mathrm{<mo\; stretchy="false">(</mo>}\lambda +\mu \mathrm{<mo\; stretchy="false">)</mo>}{u}_{j\mathrm{,}ji}+X_i\mathrm{<mo>=</mo>0,}i\mathrm{<mo>=</mo>1,2,3,}$ (2)

где u_i $—$ компоненты вектора перемещения, λ,μ $–$ коэффициенты Ламе, X_i $—$ компоненты вектора массовых сил; индексами после запятой обозначено дифференцирование по соответствующим пространственным координатам.

Пусть полупространство занимает область x₃ ≤ 0. Тогда поверхность Земли задается уравнением x₃ = 0. Эта поверхность свободна от напряжений в направлении оси x₃, следовательно на x₃ = 0 заданы граничные условия вида:

${\sigma }_{31}{\mathrm{<mo>|</mo>}}_{x_3\mathrm{<mo>=</mo>0}}\mathrm{<mo>=</mo>}{\sigma }_{32}{\mathrm{<mo>|</mo>}}_{x_3\mathrm{<mo>=</mo>0}}\mathrm{<mo>=</mo>}{\sigma }_{33}{\mathrm{<mo>|</mo>}}_{x_3\mathrm{<mo>=</mo>0}}\mathrm{<mo>=</mo>0.}$ (3)

Напряжения, создаваемые очагом готовящегося землетрясения, стремятся к нулю на бесконечности:

$\underset{x_1\to \pm \infty }{\lim }{\mathrm{\sigma }}_{\mathrm{ij}}=\underset{x_2\to \pm \infty }{\lim }{\mathrm{\sigma }}_{\mathrm{ij}}=\underset{{\mathrm{x}}_3\to -\mathrm{\infty }}{\lim }{\mathrm{\sigma }}_{\mathrm{ij}}=0$. (4)

Компоненты вектора массовых сил X, соответствующие системе, изображенной на рис. 3, выражаются следующим образом:

$X_i\mathrm{<mo>=</mo>}{p}_{ij}\frac{\partial \delta \mathrm{<mo\; stretchy="false">(</mo>}\text{x}-\xi \mathrm{<mo\; stretchy="false">)</mo>}}{\partial {\xi }_j}\mathrm{,}i\mathrm{<mo>=</mo>1,2,3,}$ (5)

где p_ij $—$ интенсивность соответствующей пары сил; δ(x−ξ) — дельта-функция; (ξ₁,ξ₂,ξ₃) — точка приложения системы сил.

Аналитическое решение. Для задачи (2) с граничными условиями (3) и (4) известны функции Грина, полученные Р. Миндлиным. Для единичной силы, приложенной к точке (ξ₁,ξ₂,ξ₃) упругого полупространства и направленной вдоль оси x₃, функция Грина g³(x) имеет вид:

$\begin{array}{l}{g}_1^3=\frac{(x_1-{\xi }_1)}{16\mathrm{\pi \mu }(1-\nu )}\left[\frac{(x_3-{\xi }_3)}{r_1^3}+\frac{(3-4\nu )(x_3-{\xi }_3)}{r_2^3}+\frac{4(1-\nu )(1-2\nu )}{r_2(r_2-x_3-{\xi }_3)}+\frac{6x_3{\xi }_3}{r_2^5}\right],\\ g_2^3=\frac{(x_2-{\xi }_2)}{16\mathrm{\pi \mu }(1-\nu )}\left[\frac{(x_3-{\xi }_3)}{r_1^3}+\frac{(3-4\nu )(x_3-{\xi }_3)}{r_2^3}+\frac{4(1-\nu )(1-2\nu )}{r_2(r_2-x_3-{\xi }_3)}+\frac{6x_3{\xi }_3(x_3+{\xi }_3)}{r_2^5}\right],\\ g_3^3=\frac{1}{16\mathrm{\pi \mu }(1-\nu )}\left[\frac{(3-4\nu )}{r_1}+\frac{5-12\nu +8{\nu }^2}{r_2}+\frac{{(x_3-{\xi }_3)}^2}{r_1^3}+\frac{(3-4\nu ){(x_3+{\xi }_3)}^2-2x_3{\xi }_3}{r_2^3}\right],\end{array}$ (6)

где ν $—$ коэффициент Пуассона, а r₁ и r₂:

$\begin{array}{l}{r}_1=\sqrt{({x_1-{\xi }_1)}^2+{(x_2-{\xi }_2)}^2+{(x_3-{\xi }_3)}^2}\mathrm{,}\\ r_2=\sqrt{({x_1-{\xi }_1)}^2+{(x_2-{\xi }_2)}^2+{(x_3+{\xi }_3)}^2}\mathrm{.}\end{array}$ (7)

Функция Грина g¹(x) для единичной силы, направленной вдоль оси x₁, выражается в виде:

$\begin{array}{l}{g}_1^1=\frac{1}{16\mathrm{\pi \mu }\left(1-\nu \right)}\left\{\frac{3-4\nu }{{\text{r}}_1}+\frac{1}{{\text{r}}_2}+\frac{{\left(x_1-{\xi }_1\right)}^2}{{\text{r}}_2^3}+\frac{\left(3-4\nu \right)\left(x_1-{\xi }_1\right)}{{\text{r}}_2^3}+\frac{4\left(1-\nu \right)\left(1-2\nu \right)\left[{\text{r}}_2^2-{\left({\text{x}}_1-{\text{ξ}}_1\right)}^2-{\text{r}}_2\left({\text{x}}_3+{\text{ξ}}_3\right)\right]}{{\text{r}}_2{\left({\text{r}}_2-{\text{x}}_3-{\text{ξ}}_3\right)}^2}\right\},\\ g_2^1=\frac{\left(x_1-{\xi }_1\right)\left(x_2-{\xi }_2\right)}{16\mathrm{\pi \mu }\left(1-\nu \right)}\left[\frac{1}{{\text{r}}_1^3}+\frac{\left(3-4\nu \right)}{{\text{r}}_2^3}-\frac{6x_3{\text{ξ}}_3}{{\text{r}}_2^5}-\frac{4\left(1-\nu \right)\left(1-2\nu \right)}{{\text{r}}_2{\left({\text{r}}_2-{\text{x}}_3-{\text{ξ}}_3\right)}^2}\right],\\ g_3^1=\frac{\left(x_1-{\xi }_1\right)}{16\mathrm{\pi \mu }\left(1-\nu \right)}\left[\frac{\left(x_3-{\xi }_3\right)}{{\text{r}}_1^3}+\frac{\left(3-4\text{ν}\right)\left({\text{x}}_3-{\text{ξ}}_3\right)}{{\text{r}}_2^3}-\frac{4\left(1-\nu \right)\left(1-2\nu \right)}{{\text{r}}_2\left({\text{r}}_2-{\text{x}}_3-{\text{ξ}}_3\right)}-\frac{6x_3{\text{ξ}}_3\left({\text{x}}_3+{\text{ξ}}_3\right)}{{\text{r}}_2^5}\right],\end{array}$ (8)

Ввиду симметричности задачи, функция Грина g²(x), соответствующая действию единичной силы вдоль оси x₂, может быть получена из функции g¹(x) заменой осей x₁ и x₂. Функции Грина, отвечающие действию двойных сил, могут быть получены дифференцированием функций gⁱ(x) по пространственным координатам, то есть в виде ∂gⁱ(x)/∂x_j [12]. При i = j получим решение для пары двойных сил, направленных вдоль соответствующей оси, при i≠ j $—$ для пары двойных сил, направленных вдоль оси i с моментом относительно оси с номером, отличным от i, j.

В общем случае, решения для смещений в упругом полупространстве можно получить при помощи формулы Вольтерра [13]:

$u_k\mathrm{<mo\; stretchy="false">(</mo>}\text{x}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_{\Sigma }}{s}_i\mathrm{<mo\; stretchy="false">(</mo>}\xi \mathrm{<mo\; stretchy="false">)</mo>}{\sigma }_{ij}^k\mathrm{<mo\; stretchy="false">(</mo>}\xi \mathrm{,}\text{x}\mathrm{<mo\; stretchy="false">)</mo>}{n}_{j}d\Sigma \mathrm{,}$ (9)

где s_i(ξ) $—$ смещение на поверхности разрыва Σ, n_j $—$ компоненты единичного вектора нормали к поверхности Σ.

Учитывая, что напряжения могут быть выражены через деформации в соответствии с законом Гука, формула Вольтерра для случая однородной и изотропной среды может быть записана в виде:

$u_k\mathrm{<mo\; stretchy="false">(</mo>}\text{x}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle {\int }_{\Sigma }}\left[\mu \mathrm{<mo\; stretchy="false">(</mo>}{s}_p{n}_q+s_q{n}_p\mathrm{<mo\; stretchy="false">)</mo>}+\lambda s_k{n}_k{\delta }_{pq}\right]\frac{\partial g_k^p\mathrm{<mo\; stretchy="false">(</mo>}\text{x}\mathrm{,}\xi \mathrm{<mo\; stretchy="false">)</mo>}}{\partial {\xi }_q}d\Sigma \mathrm{<mo>=</mo>}{\displaystyle \underset{\Sigma }{\int }}{m}_{pq}\frac{\partial g_k^p\mathrm{<mo\; stretchy="false">(</mo>}\text{x}\mathrm{,}\xi \mathrm{<mo\; stretchy="false">)</mo>}}{\partial {\xi }_q}d\Sigma \mathrm{,}$ (10)

где $m_{pq}\mathrm{<mo>=</mo>}\mu \mathrm{<mo\; stretchy="false">(</mo>}{s}_p{n}_q+s_q{n}_p\mathrm{<mo\; stretchy="false">)</mo>}+\lambda s_k{n}_k{\delta }_{pq}$ $—$ тензор плотности сейсмического момента [11], который отражает механику очага землетрясения.

Таким образом, в случае точечного источника, решение поставленной задачи может быть найдено в следующем виде:

$u_k\mathrm{<mo\; stretchy="false">(</mo>}\text{x}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{m}_{pq}\frac{\partial g_k^p\mathrm{<mo\; stretchy="false">(</mo>}\text{x},\xi \mathrm{<mo\; stretchy="false">)</mo>}}{\partial {\xi }_q}\mathrm{.}$ (11)

Формула (11) будет использована в дальнейшем как основная при моделировании.

Поскольку связь между компонентами тензора деформации ε_ij и потенциальной энергией упругих деформаций E квадратична, то повышающий коэффициент, позволяющий рассчитать напряженно-деформированное состояние земной коры при подготовке землетрясения, будем полагать равным η^−0,5. Удобный, с точки зрения вычислений, вариант оценки этого коэффициента был дан И. П. Добровольским [15]:

$\eta {\mathrm{<mo>=</mo>10}}^{\mathrm{0,26}{M}_W-\mathrm{3,93}}\mathrm{,}$ (12)

где M_W $—$ моментная магнитуда землетрясения.

Результаты моделирования

Исходные данные. Для моделирования из каталога механики очагов землетрясений «The Global Centroid-Moment-Tensor Catalog» [14] была отобрана информация обо всех землетрясениях, произошедших вблизи полуострова Камчатка с 1976 г. по 2020 г. В частности использовались данные о дате и времени землетрясения, координатах эпицентра, глубине гипоцентра, магнитуде землетрясения, тензоре плотности сейсмического момента, скалярном сейсмическом моменте. Координаты эпицентров землетрясений находились в диапазоне широт от 49° с. ш. до 60° с. ш. и диапазоне долгот от 149° в. д. до 170° в. д. Всего выборка содержала N = 877 землетрясений.

Результаты наблюдений показывают, что ориентация акустического излучения зависит от азимутального направления на эпицентры землетрясений [7]. Поэтому все рассматриваемые землетрясения были разделены на три кластера методом K-средних по пространственному расположению их эпицентров (рис. 4a).

 

Рис. 4. Все N = 877 землетрясений, разделенные на кластеры по пространственному расположению их эпицентров: •$—$ кластер № 1, •$—$ кластер № 2, • $—$ кластер № 3 (a). Землетрясения, которые по результатам моделирования могли вызвать повышенные деформации в пункте наблюдения «Микижа» $—$ ■ (b).

[Figure 4. All N = 877 earthquakes divided into clusters with respect to spatial locations of their epicenters: •$—$ cluster № 1, •$—$ cluster № 2, •$—$ cluster № 3. Earthquakes which, according to modeling results, could cause deformation intensification at Mikizha observation site $—$ ■ (b).]

 

Моделирования направления осей главных напряжений. Для каждого землетрясения из выборки было произведено моделирования напряженно-деформированного состояния земной коры, возникающее при его подготовке. Рассчитаны деформации

${\epsilon }_{max}\mathrm{<mo>=</mo>}\frac{1}{2\mu }{\sigma }_{\max }\mathrm{,}$ (13)

где ${\sigma }_{\max }\mathrm{<mo>=</mo>}\max \mathrm{<mo>\{</mo><mo>|</mo>}{\sigma }_1-{\sigma }_2\mathrm{<mo>|</mo>,<mo>|</mo>}{\sigma }_2-{\sigma }_3\mathrm{<mo>|</mo>,<mo>|</mo>}{\sigma }_1-{\sigma }_3\mathrm{<mo>|</mo><mo>\}</mo>}$ $—$ максимальное касательное напряжение, а σ₁,σ₂,σ₃ $—$ главные значения тензора напряжений.

В результате были отобраны только те землетрясения, которые вызывают повышенные предсейсмические деформационные возмущения в пункте наблюдений «Микижа» (рис. 4b). Всего таких землетрясений было 169, большая часть из них принадлежат кластеру № 2. Далее для каждого землетрясения в точке, соответствующей пункту наблюдений, вычислены собственные векторы тензора напряжений, определяющие направления главных осей. Ось максимального сжатия сонаправлена с главной осью тензора напряжений, соответствующей максимальному из значений |σ₁|, |σ₂|. Гистограммы распределения направления этой оси представлены на рис. 5. Цвет гистограммы соответствует кластеру землетрясений. Серым на гистограммах обозначена область ориентации оси максимального сжатия в пункте наблюдений «Микижа» по результатам анализа аномалий направленности геоакустической эмиссии.

 

Рис. 5. Гистограммы распределения направлений оси максимального сжатия в пункте наблюдений «Микижа» при подготовке землетрясений из кластеров № 1 (a), № 2 (b), № 3 (c) по результатам моделирования,  ■$—$ область ориентации оси максимального сжатия по результатам анализа аномалий направленности геоакустической эмиссии.

[Figure 5. Distribution histograms of main stress axis direction angles at Mikizha observation site during preparation of the earthquakes from the clusters № 1 (a), № 2 (b), № 3 (c) according to the modeling results,  ■$—$ orientation region of the compression axis according to the analysis results of geoacoustic emission directivity anomalies.]

 

Как видно из рис. 5, результаты моделирования осей главных напряжений при подготовке землетрясений из кластеров № 1 и № 2 наиболее согласуются с экспериментальными оценками. Модальный интервал гистограммы для землетрясений из кластера № 1 имеет границы от 300° до 315° и полностью покрывается интервалом экспериментальной оценки направления оси максимального сжатия (рис. 5a). Границы модального интервала для землетрясений из кластера № 2 лишь частично покрываются этим интервалом. Кроме этого присутствуют локальные максимумы значений гистограммы распределения в области от 30° до 75°, которые перпендикулярны направлению оси максимального сжатия, полученному в результате экспериментальной оценки (рис. 5b). Направления оси, полученные в результате моделирования напряженно-деформированного состояния земной коры при подготовке землетрясений из кластера № 3, плохо согласуются с экспериментальными оценками (рис. 5c). Это может быть связано с тем, что при оценке оси максимального сжатия по аномалиям направленности геоакустической эмиссии, землетрясения из этого региона не рассматривались [8].

Заключение

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

В работе использована модель, разработанная ранее в рамках линейной теории упругости [16]. Земля рассмотрена в виде однородного изотропного упругого полупространства, очаг готовящегося землетрясения $—$ в виде комбинации двойных пар сил. В модели также учтена потенциальная энергия упругих деформаций, накопленная в процессе подготовки землетрясения.

В рамках этой модели рассчитаны ориентации осей главных напряжений в пункте наблюдений «Микижа». Использованы данные из каталога механики очагов землетрясений «The Global Centroid-Moment-Tensor Catalog» о сейсмических событиях, произошедших вблизи полуострова Камчатка с 1976 по 2020 годы. Все рассматриваемые землетрясения были разделены на три группы по пространственному расположению их эпицентров.

Показано, что гистограммы распределения направлений осей главных напряжений согласуются с результатами экспериментальных оценок для двух групп землетрясений. Модальные интервалы и локальные максимумы гистограмм распределения попадают в диапазоны оценок ориентации осей максимального сжатия от 290° до 320° и от 20° до 50° соответственно.

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

Об авторах

Максим Игоревич Гапеев

Институт космофизических исследований и распространения радиоволн ДВО РАН

Автор, ответственный за переписку.
Email: gapeev.sci@yandex.ru
ORCID iD: 0000-0001-5798-7166

младший научный сотрудник лаборатории акустических исследований

Россия, 684034, c. Паратунка, ул. Мирная, д. 7

Александра Андреевна Солодчук

Институт космофизических исследований и распространения радиоволн ДВО РАН

Email: gapeev.sci@yandex.ru
ORCID iD: 0000-0002-6761-8978

кандидат физико-математических наук, старший научный сотрудник лаборатории акустических исследований, ученый секретарь

Россия, 684034, c. Паратунка, ул. Мирная, д. 7

Список литературы

Дополнительные файлы