Концепция моделирования гидродинамических процессов на шахтных полях

Концепция моделирования гидродинамических процессов на шахтных полях

РЕФЕРАТ

Концепция моделирования гидродинамических процессов на шахтных полях

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

Рис. 1. Принципы классификации моделей гидрогазодинамических (ГГД) процессов в ГТС горнопромышленного региона

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

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

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

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

Ряд параметров и расчетных характеристик задаются в виде многомерных случайных величин или случайных полей. Для таких параметров в любой точке массива P в момент времени t выделяются закономерная  и случайная x¢ составляющие

.      (1)

Пусть X(t) - случайный многомерный процесс, описывающий пространственно-временные изменения ГГД параметров или расчетных характеристик, X_i=X(t_i) - фактические значения этих параметров в моменты времени t_i, X_i,m - измеренные в системе опорных точек - пунктов наблюдения с некоторой погрешностью в те же моменты, X_i,c - значения, рассчитанные на основании моделирования. Задача оптимального прогнозирования на основе неполных исходных данных состоит в минимизации математического ожидания погрешностей прогнозируемых величин на интервале времени :

.      (2)

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

В табл. 2 приняты следующие обозначения: H - напор (уровень) подземных вод; H_m - уровень шахтных вод; Q - водоприток в шахту; Q_tr - расход перебрасываемых вод из одной шахты в другую; C_gw, C_m и C_gr - концентрации веществ в подземных, шахтных водах и в поровом растворе верхних слоев грунта; q_gw и q_sw - массопотоки в грунтовые и в поверхностные воды; P_w - давление рудничного газа в выработках; Q_g - его поток через массив и систему вентиляции; C_atm - концентрация газообразных примесей и пыли в атмосфере. Обозначения для параметров, задаваемых для расчета, следующие: K - коэффициент фильтрации; d - раскрытие трещин, a и b - азимут и угол их падения; n и V_f - пористость (трещиноватость) и объем трещин в части массива; m - мощность водоносных и водоупорных слоев; H_b - напор на границе области фильтрации; e - инфильтрация; D_f и D_g - коэффициенты диффузии водных растворов и газов в горных породах; n_e и l - параметры сорбции и нейтрализации (распада) веществ в горных породах; r_g - плотность газа в выработках; y - водонасыщенность пород, V_w - объем выработанного пространства, d_w - размер выработок; p_g - газоносность пород; q_ex - темпы добычи углесодержащих пород; j - концентрация метана в рудничном газе, C_s, q_s - концентрация веществ в шахтных породах и скорость их выщелачивания; u - скорость ветра, a_g - параметр осаждения и поглощения примесей из атмосферы на поверхности земли; d_p - гранулометрический состав пылевых выбросов; V - скорость фильтрации.

Таблица 2. Взаимосвязь основных параметров ГГД процессов

В табл. 2. также обозначено: «!» - достоверно определяемый параметр, «?» - распределенный параметр или случайная величина (функция, процесс), «+» - статистически оцениваемый параметр; «-» - определение не обязательно; СД - возможность сопоставления с исходными данными: «Н» - низкая, «С» - средняя, «В» - высокая

На этапе построения гидрогазодинамических моделей решаются следующие основные задачи: 1) разбиение массива на структурные элементы (слои, блоки), 2) построение полей проницаемости, пористости, 3) задание пустотности (пространственно-временного распределения выработок по объему), 4) задание инфильтрационных водопритоков, 5) определение граничных условий и контактов с соседними шахтными полями (в том числе переток через сбойки), 6) определение методики расчета водопритока, как на уровне шахтного поля, так и в масштабе отдельных горизонтов и блоков, 7) учет технологий отработки (особенности водоотлива и дренажа, консервации шахты), 8) учет переменной газоносности пород, различных путей газовыделения и миграции газа, 9) стыковка моделей разных масштабов, 10) стыковка гидродинамических, миграционных и газодинамических моделей на границах твердой, жидкой и газообразной сред, 11) учет взаимного влияния процессов в различных фазах.

Моделирование ГГД процессов в подработанном массиве осложняется слоистостью, наличием зон водопроводящих трещин, зон сдвижения, последовательностью отработки. Адекватный учет этих особенностей предполагает пространственную дискретизацию массива. Размер блока дискретизации должен допускать осреднение параметров и их оценку (измерение) в натурных условиях. С этой точки зрения блоки размером менее 50 м по горизонтали и 15-20 м по вертикали нецелесообразны: при этом размеры пустот (выработок) становятся соизмеримыми с размерами блока. Напротив, увеличение блоков, соответствующее региональным моделям, в большей мере соответствует гипотезе сплошности среды, на которой строится модель фильтрации. При этом напор, рассчитанный на крупномасштабных моделях, особенно в массиве с крутопадающими слоями при наличии осложняющих факторов, следует трактовать как средне вероятный в пределах блока.

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

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

а)                                                               б)

Рис. 2. Изменение коэффициента фильтрации K (м/сут) с глубиной z на полях шахт: а) Александр-Запад, б) им. Калинина; крайние значения; - средние значения; N - число измерений (столбцы)

Анализ графиков свидетельствует о большой вариации K в пределах одного слоя, что характерно для всех глубин. Это может быть связано: 1) с наличием нескольких подзон различной проницаемости, 2) с макронеоднородностью и масштабными эффектами, 3) с хаотической неоднородностью. Значительный разброс данных измерений подтверждает необходимость использования моделей, в которых проницаемость пород является случайным полем или распределенным параметром.

Результаты измерений, представленные на рис. 2, можно рассматривать как косвенную оценку надежности исходных данных для моделирования, которая возрастает с увеличением N. В то же время адекватная оценка статистической значимости K должна базироваться на полной информации о результатах опытов, в соответствии с которыми проверяется гипотеза о виде распределения и его моментах.

Закономерная неоднородность на уровне отдельных блоков может рассматриваться как стохастическая в пределах слоя (горизонта) или его больших подзон. В разработанных моделях задается пространственно-временная изменчивость параметров, связанная с неодновременным проведением горных работ на различных горизонтах отработки.

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

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

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

Изменение проницаемости подработанного массива учитывается с помощью переменных во времени характеристик, отражающих динамику сопутствующих фильтрации процессов: просадку и уплотнение пород, суффозию, вытеснение рудничного газа и др.

Стыковка моделей миграции в грунтах и подземных водах осуществляется через параметры массопотока на границах «подземные - поверхностные воды» и на границе зоны аэрации. Модели миграции в верхних слоях грунтов и отвалах шахтных пород и воздушной среде сочленяются посредством параметров осаждения твердых частиц и поглощения газообразных примесей на подстилающей поверхности. Принципы моделирования движения пылевой и газовой компонент в геотехнических системах изложены далее в главе 3.

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

Объем пустот в пределах шахтного поля V состоит из объема трещин и пор V_f, а также объема выработанного пространства V_w. Объем V_f учитывает как естественную, так и техногенную трещиноватость. Изменение параметра V_w в пространстве задается на основе фактических планов горных работ с учетом коэффициента заполнения [3].

С целью адекватного отображения фильтрационной неоднородности и пространственного распределения зон отработки шахтного поля производится его разбивка на блоки в виде прямоугольников. На каждом горизонте толщиной Dz, где ведется отработка, задается послойно общий объем выработанного пространства V_w,l (l=1,… N_z, N_z - число слоев). При схематизации структуры массива на каждом горизонте отработки выделяется несколько зон с отличающимися значениями фильтрационных параметров и объемами выработанного пространства. Тогда, в соответствии с горизонтальным распределением объема V_w по слоям и с учетом последовательности отработки для каждого блока можно задать объем выработок, время их появления и погашения. Исходя из объема V_w,l определяется площадь горизонтального сечения выработок S_h на заданной глубине z как отношение объема пустот в слое к его толщине:  [4].

Пример такой фильтрационной схематизации для шахтного поля в условиях Центрального района Донбасса показан на рис. 3 и 4. Наличие зон усиленной инфильтрации (отстойников шахтных вод, понижений рельефа) моделируется путем задания неравномерной по площади интенсивности инфильтрационного питания e.

шахта гидродинамический массив

Рис. 3. Распределение выработанного пространства V_w и среднего коэффициента фильтрации K_f по глубине

В условиях нарушения естественных геологических структур при разработке пластов, наличия зон обрушения с повышенной вертикальной проницаемостью расчет изменения уровня подземных вод целесообразно вести на основе уравнения плановой неустановившейся фильтрации [5,6]

. (3)

Здесь Н - уровень подземных вод, T_x и T_y - проводимость вдоль осей Ox и Oy соответственно, n_f - трещинная пористость пород массива, Q_w - интенсивность стоков в выработки, e - интенсивность инфильтрации.

В сложных гидрогеологических условиях фильтрационный расчет выполняется с помощью численных моделей. Для решения уравнения (3) использована попеременно-треугольная схема метода конечных разностей, сочетающая преимущества явной схемы с хорошей вычислительной устойчивостью [7,8].

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

Рис. 4. Расположение горизонтов отработки и очагов инфильтрации в пределах шахтного поля:1,…, 5 - контуры горизонтов отработки на отметках 1 - от 75 до 210 м; 2 - от -100 до 75 м; 3 - от -300 до -100 м; 4 - от -500 до -300 м; 5 - от -600 до -500 м; S - шахтный ствол

Параметры водопроводимости определяются по формулам

, (4)

где K_x и K_y - значения коэффициента фильтрации вдоль осей Ox и Oy соответственно. При послойной аппроксимации параметра проницаемости интегрирование в выражениях (4) заменяется соответствующей суммой.

Особенностью разработанной методики является способ расчета притока в выработки в соответствии с принятой схематизацией подработанного массива. Водоприток в выработки в пределах блока (i, j) определяется суммированием водопритоков по всем дренируемым горизонтам отработки, число которых равно N_ij

.  (5)

Величина N_ij изменяется по мере разработки и дренирования более глубоких горизонтов и в зависимости от положения уровня воды в шахтном стволе, гидравлически связанном с выработками.

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

,         (6)

где K_w,ijk - коэффициент фильтрации на данном горизонте в блоке (i, j), S_w,ijk - площадь поверхности «выработки - окружающий массив», DH_w,ijk - перепад напора в водонасыщенных породах вокруг выработок и поверхностью стока, L_a,ijk - средняя длина пути фильтрации из произвольной точки массива к поверхности выработок. Параметры, входящие в выражение (6), можно оценить на основе статистических характеристик нарушенного горного массива.

Так, площадь фильтрации S_w (далее индексы i, j, k для простоты опущены) определяется следующим образом. Исходя из приближения к цилиндрической форме выработок, можно записать для объема пустот V_w в пределах некоторого блока-параллелепипеда

, (7)

где l_w(r) - длина выработки радиуса r,  - ее среднее значение для выработок всех размеров, p_w(r) - плотность распределения выработок по объему блока, [r_min, r_max] - диапазон изменения радиусов выработок.

При аналогичных допущениях

, (8)

где p_s(r) - плотность распределения выработок по площади их боковой поверхности в зависимости от радиуса. Сопоставляя (2.7) и (2.8), получим

. (9)

Параметры распределений p_w(r) и p_s(r) можно оценить при наличии информации о размерах выработок в данном блоке или части шахтного поля. При равномерном распределении (p_w(r) = p_s(r) =1/(r_max - r_min)₎

.     (10)

Тогда при обычном диапазоне радиусов выработок (r_min=0,5 м, r_max=2,2 м) S_w»1,32 V_w. С увеличением доли малых выработок числовой коэффициент перед V_w возрастает, хотя обычно не превышает 2-2,5. Увеличение площади контакта «выработка - окружающий массив» аналогично увеличению удельной поверхности скелета пористой среды при уменьшении размеров частиц [8].

При осреднении напора по вертикали в принятой расчетной схеме перепад напора DH_w,ijk рассчитывается как разность между уровнем подземных вод и средним положением дренируемых (или затапливаемых) выработок на горизонте отработки:

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

Параметр L_a в формуле (6) можно оценить следующим способом. Предположим, что пустотность (доля объема выработок в массиве) в горизонтальном или вертикальном сечении блока равна его объемной пустотности h_a= V_w/V_b, где V_b - объем блока. Тогда суммарная площадь сечения выработок через некоторую грань площадью S¢ составляет S_w,S = h S¢. Средняя площадь сечения отдельной выработки гранью блока

, (11)

где a - угол между выработкой и гранью блока. При равномерном распределении p_s(r) и тех же значениях r_max и r_min величина  составит около 9 м². Предположим далее, что отношение площади непроницаемой фазы к площади выработок в сечении остается одинаковым в масштабе блока и для более мелкого блока вокруг отдельной выработки:

,     (12)

где S₁ - площадь непроницаемой фазы в сечении, приходящейся на одну выработку.

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

.  (13)

Рис. 5. Схема к расчету водопритока в выработки из массива

При типичных значениях пустотности и размеров выработок значение r_a составляет несколько метров, уменьшаясь в сильно нарушенных участках и возрастая в частях шахтного поля с малым объемом выработанного пространства. С учетом извилистости пустот следует принять L_a=c r_a, где c - коэффициент извилистости, принимаемый по аналогии с пористой средой [9].

Изменение уровня воды в стволе в период осушения и добычи определяется по графику отработки угольных пластов. После отключения водоотлива на каждом шаге по времени величина H_s рассчитывается по формуле [4]:

 (14)

где Q_w,S(t) - суммарный приток в выработки в момент времени t, S_h(H_s) - площадь горизонтальной поверхности выработок на уровне H_s.

Разработанная методика применена для расчета изменения уровня подземных вод на шахтном поле в период его разработки, и после отключения водоотлива. Исходные данные принимались по горно-гидрогеологическим условиям, типичным для Центрального района Донбасса. На границе шахтного поля размером 5´6 км установился стационарный уровень подземных вод 200 м. Трещинная пористость пород массива составляет 0,03. При затоплении, для учета влагосодержания в трещинах, вместо этой величины в уравнении (3) задавался меньший параметр водоотдачи. Параметры, характеризующие распределение коэффициента фильтрации и объема выработанного пространства по глубине, приведены на рис. 3 и 4. Расчеты проводились на прямоугольной сетке с размером блоков 50´50 м, с максимальным шагом по времени 10 сут. Результаты моделирования представлены на рис. 6 - 8.

Анализ изменения суммарного водопритока в шахту (рис. 6) показывает, что наибольший приток наблюдается при отработке верхних слоев, характеризующихся большей проницаемостью. Этим объясняются, в частности, локальные максимумы величины Q_w в первые 10-20 лет после начала разработки. В этот период уровень подземных вод находится преимущественно в пределах двух более проницаемых верхних горизонтов отработки. По мере заглубления выработок и дренирования нижних горизонтов водоприток стабилизируется.

Рис. 6. Изменение водопритока (1), уровня воды в шахтном стволе (2) и уровня подземных вод вблизи шахтного ствола (3)

Рис. 7. Уровень подземных вод на шахтном поле перед отключением водоотлива

а)

б)

Рис. 8. Изменение уровня подземных вод в шахтном поле после отключения водоотлива: а) вдоль профиля А-А, б) вдоль профиля В-В (см. рис. 7); 1, 2, 3 - моменты времени.

Период незначительного уменьшения водопритока (28-55 лет после начала отработки) характеризуется практически постоянным уровнем воды в шахтном стволе и сформировавшейся воронкой депрессии подземных вод (рис. 7), которая повторяет контуры зон отработки. Наибольшее понижение приурочено к тем зонам, где, в соответствии с принятыми исходными данными, объем выработанного пространства больше. По мере подъема уровня подземных вод при затоплении отмеченные неравномерности сглаживаются (рис. 8).

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

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

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

1. Использование и охрана подземных вод / Н.А. Маринов, А.Е. Орадовская, Е.В. Пиннекер и др. Новосибирск, Наука, 1983. 231 с.

2. Гонтаревский В.П., Кулешов В.М. Обеспечение эколого-гидробезопасности при ликвидации шахт // Уголь Украины. - 1999. - №11. - С. 22-25.

3. ДБН Ф. 2.2-1-95. Состав и содержание материалов оценки воздействий на окружающую среду при проектировании и строительстве предприятий, зданий и сооружений. Основные положения проектирования / Госкомградостроительства Украины, Минэкобезопасности Украины. - К.: Укрархстройинформ, - 1996. - 14 с.

4. КД 12.12.004-98 «Ликвидация угольных шахт. Защита земной поверхности от затопления горных выработок. Рекомендации» / УкрНИМИ, 1998. - 46 с.

5. Инструкции по безопасному ведению горных работ у затопленных выработок / ВНИМИ, 1984. - 66 с.

6. Бошенятов Е.В. Влияние метаморфизма горных пород на величину безопасной глубины разработки под водными объектами // Горная геомеханика и маркшейдерское дело: Сб. науч тр. ВНИМИ. - СПб, 1999. - С. 150-153.

7. Панасенко Г.П. О методике определения мощности зоны водопроводящих трещиноватых пород над очистными выработками угольных шахт Восточного Донбасса // В сб. «Результаты геологических исследований на территории Нижнего Дона и Нижней Волги». - Ростов-на-Дону: Изд-во РГУ, 1967. - С. 54-61.

8. Мохов А.В. Определение высоты зоны водопроводящих трещин на каменноугольных месторождениях в натурных условиях. - М., 1988. - 25 с.

9. Беседа Н.И. Особенности гидрогеологиченских условий вскрытия и разработки углей на больших глубинах в юго-западной части Донбасса / УкрНИИНТИ. - К., 1970. -49 с.

10.   Методические указания по изучению и прогнозированию гидрогеологических условий при разработке угольных месторождений Донецкого бассейна. Днепропетровск: ДО ИМР, 1979. - 100 с.

11.   Беседа Н.И., Сляднев В.А., Яковлев Е.А. Региональные техногенные изменения геологической среды Донбасса под влиянием горных работ. - К. - 1977. - 66 с.

12.   Цабут И.И., Улицкий О.А., Шворак Н.В. Анализ многолетней динамики обводнения угольных шахт Центрального углепромышленного района Донбасса // Труды ДонГТУ, Серия горно-геологическая, вып. 11, Донецк, ДонГТУ, 2000. - С. 124-128.

13.   Ермаков В.Н., Улицкий О.А., СпожакинА.И. Изменение гидродинамического режима шахт при затоплении // Уголь Украины. - 1998. - №6. - С. 11-13.

14.   Ермаков В.Н., Улицкий О.А., Питаленко Е.И. и др. Изменение гидрогеологических и геомеханических условий при закрытии шахт // Наук. Праці ДонНТУ, Серія гірничо-геологічна, вип. 321, Донецьк, ДонНТУ, 2001. - С. 69-73.

15.   Ермаков В.Н., Семенов А.Н., Улицкий О.А. и др. Развитие процессов подтопления земной поверхности под влиянием закрывающихся шахт // Уголь Украины. - 2001. - №6. - С. 12-14

16.   ДНАОП 1.1.30 - 1.01.96. Правила безпеки у вугільних шахтах. - Донецьк: Донеччина, 2001. - 495 с.

17.   Ермаков В.Н., Петренко С.Я., Касимов О.И., Кочерга В.Н. О предотвращении выделения газов из ликвидируемых шахт Стахановского региона // Уголь Украины, - №5. - 1999. - С. 15-17.

18.   Василянский Н.П., Кочерга В.Н. Предотвращение газовыделения на земную поверхность путем дегазации выработанных пространств / Сб. науч. трудов МакНИИ. - 1987. - С. 27-32.

19.   Айкер Х., Хесбрюгге С. Газовыделение из закрытых шахт и погашенных выработок // Глюкауф. - №23. - 1984. - С. 29-34.

20.   Robinson R. Mine gas hazards in the surface environment // Mining Technology, Section A. - Vol. 109. - 2000. - P. A228-236.