Моделирование в PLAXIS. Причины оседания поверхности земли на орошаемых территориях вблизи угольного карьера

ВВЕДЕНИЕ

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

Тзапоглу и Лупасакис изложили полученные результаты в статье «Численное моделирование факторов, вызывающих оседание поверхности земли из-за чрезмерной эксплуатации водоносного горизонта на территории расположения угольного карьера Аминдео, Греция» [1], которая была опубликована в журнале HydroResearch в 2019 году. Здесь мы кратко рассмотрим ее материалы.

Чтобы подчеркнуть важность проблемы, авторы работы [1] напомнили, что чрезмерная эксплуатация водоносных горизонтов в конце XX века при удовлетворении растущих промышленных, сельскохозяйственных и бытовых потребностей в воде привела к обширным оседаниям поверхности земли более чем в 150 густонаселенных регионах мира, например в китайских Пекине, Шанхае, Тяньцзине, Сучжоу, Чанчжоу, в американских Лас-Вегасе, Неваде, Хьюстоне, Аризоне, в индонезийских Джакарте, Бандунге, в испанской Муртии, в итальянской Венеции, на греческих территориях региона Калохори, долины Антемонта, восточной и западной равнин Лариса, равнины Триасион, угольного суббассейна Аминдео и в других местах.

Тзапоглу и Лупасакис [1] остановились на исследовании медленного оседания поверхности земли, которое в последние десятилетия постоянно происходит на территориях сельхозугодий и населенных пунктов вокруг угольного карьера Аминдео (Amyntaio, греч. Αμύνταιο) в периферийной единице Флорина периферии Западная Македония на севере Греции (рис. 1). Этот процесс вызван чрезмерной откачкой грунтовых вод через дренажные скважины для осушения вышеназванного угольного карьера и через скважины системы орошения сельскохозяйственных земель на окружающей равнине. В результате оседание распространилось уже на расстояние 4 км вокруг карьера и нанесло значительный ущерб населенным пунктам Анаргири (Anargiroi), Валтонера (Valtonera), Фанос (Fanos) и другим, а также дорожной сети.

Рис. 1. Спутниковый снимок угольного карьера Аминдео (Amyntaio Mine) и прилегающей к нему территории [1]. Черным шрифтом написано название карьера, белым шрифтом – названия населенных пунктов, синим – названия озер

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

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

Основные цели работы [1] – исследование механизма оседания поверхности земли в рассматриваемом районе и оценка степени влияния на него основных факторов. 

Для достижения этих целей авторы статьи [1] использовали не только анализ данных из литературных источников, полевые и лабораторные измерения, но и (прежде всего) численное моделирование в программном комплексе PLAXIS 3D на основе метода конечных элементов.  

ИСХОДНЫЕ ДАННЫЕ ДЛЯ МОДЕЛИРОВАНИЯ

Все необходимые для моделирования исходные данные, касающиеся инженерно-геологических условий, Тзапоглу и Лупасакис [1] собрали на основе анализа результатов более ранних исследований рассматриваемой территории (собственных и других авторов), а также новых собственных полевых измерений и наблюдений. В том числе они выполнили картирование разрывов поверхности земли на исследуемой территории общей протяженностью 15,5 км (с помощью системы глобального позиционирования), а также в мае 2016 года провели полевые измерения уровней грунтовых вод (УГВ) в водоносном горизонте с использованием имеющихся скважин, пробуренных ранее Институтом геологии и эксплуатации недр (IGME), и нескольких новых скважин.

Район исследований относится к пелагонской (пелагонийской) геотектонической зоне Греции и сложен альпийско-предальпийскими, неогеновыми и четвертичными геологическими формациями (рис. 2, табл. 1).

Рис. 2. Геологическая карта угленосного суббассейна Аминдео (Аминтайо, Amyntaio) [1]. См. также таблицу 1

Таблица 1. К условным обозначениям на рисунке 2. – Ред.

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

Неогеновые геологические формации подразделяются на нижнюю и верхнюю свиты, то есть на лигнитсодержащую комнинскую (Komnina + Vegora) свиту, которую слагают глины и пески с прослоями песчаников, и лигнитсодержащую птолемаидскую (Ptolemaida) свиту, сложенную в основном переслаивающимися пылеватыми глинами и глинистыми мергелями (но с присутствием и слоев лигнита). Эти неогеновые формации прочные, переуплотненные и водонепроницаемые. На них не могут повлиять никакие процессы, вызывающие дальнейшее уплотнение. Кровли их верхних слоев представляют собой подошвы водоносных горизонтов.

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

Проастиосская (Proastio) свита сложена речными отложениями нижне-среднеплейстоценового возраста (залегающими на лигнитсодержащих неогеновых формациях, а в восточной части суббассейна – под пердикасской свитой). Это переслаивающиеся гравийно-песчаные отложения (sandy gravels), красные глины, пылевато-глинистые пески (silty clayey sands) и пылеватые глины с большим количеством органики (organic silty clays). Верхние горизонты содержат более 80% гравия с песком, а глубже увеличивается до 30% содержание слоев песчанисто-пылеватых глин (sandy silty clays). Среднее значение индекса компрессии С_с (compression index) достигает 0,16 для пылевато-глинистого песка (silty clayey sand), 0,18 для пылеватой глины с большим количеством органики (organic silty clay) и 0,16 для песчанисто-пылеватой глины (sandy silty clay). Таким образом, с учетом напряжений предварительного уплотнения и в целом малого присутствия мелкозернистых слоев авторы работы [1] характеризуют отложения этой свиты как нормально уплотненные или немного переуплотненные и считают, что их вклад в общее оседание поверхности земли минимален.

Пердикасская (Perdikas) свита залегает на большей части исследуемой территории под молодыми аллювиальными отложениями. Сложена озерными и речными отложениями нижне-среднеплейстоценового возраста и состоит из перемежающихся слоев глинистых грунтов (песчанистой глины, глинистого песка – sandy clay, clayey sand) и песка с включениями мергеля (sand with marls), а в верхних слоях встречаются прослои пылеватых глин с большим количеством органики (organic silty clays). Крупнозернистые горизонты по мощности составляют только 10–20% стратиграфического разреза, слои песчанисто-пылеватых глин (sandy silty clays) – 65–75%, а слои пылеватых глин с большим количеством органики (organic silty clays) – 10–20%. Средние значения индекса компрессии С_с для этой свиты высоки и достигают 0,99 для пылеватой глины с большим количеством органики (organic silty clay), 0,27 для песчанисто-пылеватой глины (sandy silty clay) и 0,164 для пылевато-глинистого песка (silty clayey sand). Поэтому авторы работы [1] делают вывод, что из-за своих геотехнических характеристик и давления предуплотнения отложения этой свиты являются нормально уплотненными и подвержены осадкам при чрезмерной эксплуатации водоносного горизонта.

Анаргирийская (Anargiroi) свита сложена флювиальными (fluviotorrential) отложениями среднеплейстоценового возраста. Она состоит из перемежающихся гравийно-песчаных отложений (sandy gravels), слоев красных глин и пылевато-глинистых песков (silty clayey sands), но в более глубоких горизонтах преобладают прослои песчанисто-пылеватых глин (sandy silty clays). Крупнозернистые горизонты составляют по мощности 82–88% свиты, а глинистые прослои – от 12 до 18%. Среднее значение индекса компрессии С_с является низким, достигая 0,16 для песчанисто-пылеватой глины (sandy silty clay). Поэтому, как считают авторы статьи [1], даже учитывая расположение данной свиты над кровлей водоносного горизонта, ясно, что ее присутствие не может влиять на оседание поверхности земли.

Молодые отложения включают старые каменистые осыпи, осыпные конусы выноса, травертины (известковые туфы) и аллювиальные конусы выноса верхнеплейстоценового возраста. Кроме того, центральную часть равнины занимают тонкие слои торфа и молодые озерные отложения. Завершают стратиграфический разрез элювиальные и современные аллювиальные отложения. Индексы компрессии С_с составляют 0,187 для озерных отложений и 1,707 для торфяных. Авторов работы [1] указывают, что эти отложения, хотя и сильно подвержены уплотнению, но могут давать осадки только в результате окисления.

Геотектонические условия на севере Греции контролируются полем деформаций растяжения, активным начиная со среднего-верхнего миоцена. Разломы на территории угленосного суббассейна Аминдео можно разделить на две группы, затрагивающие в том числе четвертичные геологические формации – с направлениями от СЗ–ЮВ до ССЗ–ЮЮВ и с направлениями от СВ–ЮЗ до ВСВ–ЗЮЗ (см. рис. 2). Район исследований практически окружен активными тектоническими разломами. Три из них показаны линиями F_pxa, F_ang и F_veg на рисунке 2 (похоже, авторы рисунка ошиблись с условными обозначениями на рисунке 2 или с описанием, поскольку, судя по карте, разлом F_ang не является активным, что противоречит тексту. – Ред.).

Что касается гидрогеологических условий исследуемой территории, то серия водоносных горизонтов (которую авторы работы [1] чаще называют просто водоносным горизонтом, что в дальнейшем будем делать и мы в этом обзоре) находится в пределах разреза четвертичных отложений на глубине около 500–600 м, имеет мощность в центре угленосного суббассейна Аминдео 120–130 м и пополняется за счет притока из окружающих альпийских и предальпийских геологических формаций и карстовых водоносных горизонтов, а также за счет инфильтрации поверхностных вод.

Тзапоглу и Лупасакис [1] с помощью программ Surfer и Arc GIS 10.1 построили пьезометрическую карту рассматриваемой территории по состоянию на май 2016 года (рис. 3) (и уже имелась аналогичная карта для мая 1992 года). На основе сравнения пьезометрических уровней для мая 1992 года и для мая 2016 года на карту исследуемого района были нанесены изолинии снижения УГВ (рис. 4). Из рисунка 4 хорошо видно, что при приближении к периметру угольного карьера разница в пьезометрических уровнях за 24 года в целом становится больше, достигая максимальных значений от 70 до 100 м. Можно также заметить, что явное влияние карьера на снижение уровней УГВ распространяется на 4 км и более от его периметра и достигает берега озера Химадитис (Химадион, Химадитида – Cheimaditida) и населенных пунктов Анаргири (Anargiroi), Валтонера (Valtonera), Фанос (Fanos) и других.

Рис. 3. Изолинии пьезометрических уровней (м) по состоянию на май 2016 года, нанесенные на спутниковый снимок исследуемой территории желтым цветом [1]

Рис. 4. Изолинии снижения пьезометрических уровней (м) с мая 1992 года по май 2016 года, нанесенные на спутниковый снимок исследуемой территории желтым цветом [1]

Результаты измерений пьезометрических уровней в точках P1 и P2 (см. рис. 4) с 21 сентября 1982 года по 21 сентября 2014 года, представленные на рисунке 5 в графической форме, показывают, что для более близкой к угольному карьеру точки P1 хорошо заметна тенденция к снижению уровня грунтовых вод за 32 года. Более подробные литературные данные и полевые измерения показали, что уровень воды в водоносном горизонте на территории угленосного суббассейна Аминдео в последние десятилетия действительно имеет устойчивую тенденцию к снижению (скорее всего связанную с непрерывным потоком грунтовых вод в сторону карьера из-за работы скважин, осуществляющих его дренаж), но имеет сезонные колебания (вероятнее всего связанные с чрезмерной откачкой воды из скважин оросительной системы в сезон полива сельскохозяйственных культур).

Рис. 5. Результаты измерений пьезометрических уровней (м) с 21 сентября 1982 г. по 21 сентября 2014 г. в точках P1 и P2 (показанных красным цветом на рисунке 2) [1]. Для более близкой к угольному карьеру точки P1 хорошо заметна тенденция к снижению уровня грунтовых вод за 32 года

Авторы работы [1] приводят очень удачное сравнение: в настоящее время угольный карьер Аминдео работает наподобие гигантской скважины на воду, в результате чего образуется депрессионная воронка, простирающаяся на километры вокруг. Последняя также включает в себя более мелкие динамические депрессионные воронки, образовавшиеся при работе скважин для оросительных систем. Например, на пьезометрической карте по состоянию на май 2016 года (рис. 4) лучше всего просматриваются две такие воронки – к ВЮВ от Валтонеры (Valtonera) и к ЮЮВ от Фаноса (Fanos).

ПОСЛЕДСТВИЯ ЧРЕЗМЕРНОЙ ЭКСПЛУАТАЦИИ ВОДОНОСНОГО ГОРИЗОНТА

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

На основе подробной карты, составленной авторами статьи [1] по результатам обширных полевых исследований, они разделили разрывы поверхности на две категории (рис. 6). Большинство из них имеют направление СВ–ЮЗ и параллельны тектоническим разломам, затрагивающим четвертичные геологические формации. Остальные имеют направление СЗ–ЮВ и параллельны тектоническим разломам, затрагивающим предальпийские геологические формации. На этой основе Тзапоглу и Лупасакис [1] делают вывод, что совершенно очевидна корреляция разрывов поверхности с пространственным распределением тектонических разломов.

Рис. 6. Пространственное распределение на территории угленосного суббассейна Аминдео разрывов поверхности [1]. Черные штриховые линии на космическом снимке показывают возможные разрывы поверхности, а черные сплошные линии – точно установленные разрывы поверхности. На фотографиях вокруг спутникового снимка показаны повреждения домов и объектов инфраструктуры в местах разрывов, отмеченных на снимке красными точками и латинскими буквами A, B, C, D, E, F, G

Примерно 15,5 км зарегистрированных разрывов поверхности нанесли значительный ущерб жилым домам и объектам инфраструктуры, в том числе дорожной сети. При этом в основном пострадали населенные пункты Анаргири, Валтонера и Фанос.

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

Один основной разрыв поверхности с направлением СВ–ЮЗ пересек всю территорию населенного пункта Валтонера (который находится примерно в 1,5 км от карьера), вызвав серьезные повреждения зданий и дорог.

В Фаносе, расположенном на достаточно большом расстоянии от карьера, выявлен один разрыв поверхности с направлением СВ–ЮЗ, нанесший ущерб селению на всем своем протяжении, повредив несколько зданий.

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

Самые большие вертикальные смещения по разрывам поверхности достигают 0,5–1 м и наблюдаются близко к угольному карьеру. При удалении от него такие деформации уменьшаются. На основе этих наблюдений и карты снижения уровней воды в водоносном горизонте (см. рис. 4) авторы работы [1] делают вывод, что осадки поверхности земли в целом пропорциональны расстоянию от карьера и, следовательно, снижению УГВ.

МОДЕЛИРОВАНИЕ В ПРОГРАММНОМ КОМПЛЕКСЕ PLAXIS

Для устранения неясностей в отношении механизма оседания поверхности земли и выявления факторов, особенно сильно влияющих на его активизацию и к тому же не одинаковых в разных зонах исследуемой территории, требовались дальнейшие исследования. Тзапоглу и Лупасакис [1] выполнили серию параметрических расчетов с использованием программного комплекса PLAXIS 3D (версии 2016 года) на основе метода конечных элементов (МКЭ). Целью была оценка влияния на активизацию осадок вокруг угольного карьера таких факторов, как снижение уровня водоносного горизонта, инженерно-геологические характеристики геологических формаций, амплитуды смещений по разломам (fault offsets – скорее всего имеются в виду вертикальные амплитуды, но авторы работы [1] это не уточняют. – Ред.) и углы падения плоскостей разломов (fault angles). Исходные параметры для введения в программу (коэффициент Пуассона v; удельный вес неводонасыщенного грунта γ_unsat, кН/м³; удельный вес водонасыщенного грунта γ_sat, кН/м³; удельное сцепление C, кН/м²; угол внутреннего трения φ, град.; коэффициенты фильтрации по направлениям трех осей координат K_x, K_y, K_z, м/сут; одометрический модуль деформации E_oed, кН/м²; контрольный модуль жесткости при контрольном всестороннем давлении E_ref, кН/м²) были получены в результате оценки результатов предыдущих исследований авторов статьи [1].

При этом для моделирования влияния снижения уровня воды в водоносном горизонте на величины осадок Тзапоглу и Лупасакис [1] сделали определенные допущения:

• пердикасская (Perdikas) и проастиосская (Proastio) свиты были разделены на верхнюю и нижнюю части (из-за того, что значения индекса компрессии С_с увеличиваются с глубиной);
• было принято, что крупнозернистые и мелкозернистые слои пердикасской (Perdikas) и проастиосской (Proastio) свит имеют равный вклад в стратиграфию;
• расчеты выполнялись в шесть этапов при последовательном снижении уровня грунтовых вод на величины 10, 20, 30, 40, 50 и 60 м;
• моделирование проводилось в режиме расчета консолидации (осадки во времени) с использованием идеально-упругопластической модели Мора – Кулона.

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

Моделирование для Анаргири

Геотехнические данные для слоев грунта, использованные при моделировании для Анаргири, были получены при лабораторных испытаниях многочисленных образцов, собранных в ходе предыдущих исследований авторов статьи [1], а также путем оценки результатов, опубликованных другими учеными.

Взятый для построения модели участок территории Анаргири, расположенный в 600 м от периметра карьера, пересекается разломом F_ang (см. рис. 2), имеющим амплитуду смещения 50 м. Как видно из рисунка 7, на северо-западной стороне (справа от разлома) между анаргирийской (Anargiroi) и проастиосской (Proastio) свитами залегает пердикасская (Perdikas) свита, тогда как на юго-восточной стороне (слева от разлома) последняя отсутствует.

Рис. 7. Упрощенная имитационная 3D-модель грунтовых условий для рассмотренного участка в селении Анаргири [1]. Мощности слоев указаны в метрах. Геотехнические параметры слоев приведены в таблице 2

Таблица 2. Параметры грунтовых слоев пердикасской (Perdikas), анаргирийской (Anargiroi) и проастиосской (Proastio) свит на изучаемом участке населенного пункта Анаргири [1]

Результаты численного моделирования вертикальных смещений в грунтовом массиве из-за снижения уровня воды в водоносном горизонте в Анаргири представлены на рисунке 8. Эти деформации увеличиваются пропорционально падению УГВ. Видно, что осадки максимальны в северо-западной (правой) части Анаргири, где в разрезе имеется пердикасская (Perdikas) свита, на которую больше всего влияет снижение уровня грунтовых вод.

Рис. 8. Результаты моделирования вертикальных смещений в грунтовом массиве из-за снижения уровня грунтовых вод на рассмотренном участке населенного пункта Анаргири [1]. Цветовая шкала показывает величины вертикальных смещений (м). В таблице в углу рисунка в левом столбце приведены поэтапные величины снижения УГВ (м), а в правом столбце – соответствующие вертикальные смещения точки B относительно точки A (м), отмеченных на модели черным цветом

При поэтапном понижении уровня грунтовых вод на 10, 20, 30, 40, 50 и 60 м сравнивались величины оседания поверхности земли в точке B относительно точки А (см. рис. 8). Из таблицы в левом углу рисунка 8 видно, что оседание увеличивается со снижением уровня воды в водоносном горизонте. Наибольшая осадка в этом случае составила 0,841 м при падении УГВ на 60 м. Это довольно высокое значение, и оно говорит о том, что снижение УГВ затрагивает в основном только одну сторону от разлома (правую).

Для подтверждения достоверности результатов моделирования, полученных в PLAXIS 3D, авторами работы [1] были выполнены полевые измерения вертикальных смещений по разрывам поверхности в Анаргири (на более широкой территории, чем участок моделирования), которые оказались хорошо согласованными с итогами анализа методом конечных элементов. Некоторые из результатов измерений представлены на рисунке 9. Видно, что чем больше снижение УГВ под той или иной точкой (и чем она ближе к периметру угольного карьера), тем больше вертикальное смещение в месте разрыва (40–50, 35–45, 15–25, 10–15 см в точках A, D, C, B соответственно). Сами измеренные вертикальные смещения по разрывам были близки к полученным при моделировании. 

Рис. 9. Увеличенный космический снимок территории Анаргири, на который нанесены линии разрывов поверхности (желтым цветом) и изолинии снижения УГВ с мая 1992 года по май 2016 года (белым цветом) [1]. На фотографиях вокруг снимка желтым цветом показаны величины вертикальных смещений (см) в точках разрывов A, B, C, D, обозначенных на снимке красным цветом

Чтобы исследовать влияние основных причинных факторов на величины осадок поверхности земли, Тзапоглу и Лупасакис [1] провели ряд параметрических расчетов для трех сценариев моделирования.

При первом сценарии при прочих равных условиях оценивалось влияние соотношений мощностей сжимаемых отложений пердикасской (Perdikas) и проастиосской (Proastio) свит. Было выполнено 5 серий расчетов для следующих процентных соотношений мощностей слоев «песок / пылеватая глина»: 100/0; 70/30; 50/50; 30/70; 0/100. При этом для каждого из этих пяти случаев вводилась амплитуда смещения по разлому 50 м, угол падения плоскости разлома 65 град. и поэтапно увеличивалось снижение уровня воды в водоносном горизонте – на 10, 20, 30, 40, 50 и 60 м. На рисунке 10, А в графическом виде показаны результаты расчетов – кривые А1, А2, А3, А4, А5, соответствующие вышеуказанным соотношениям «песок / пылеватая глина» (в том же порядке). По оси абсцисс отложены величины снижения уровня грунтовых вод (м), а по оси ординат – оседание поверхности земли (м) в точке В относительно точки А (см. рис. 8) со знаком «минус». Из рисунка 10, А видно, что чем больше содержание пылевато-глинистых слоев, тем ниже идет кривая, то есть тем больше оседание поверхности земли.

При втором сценарии параметрического моделирования при прочих равных условиях исследовалось влияние амплитуды смещения по разлому. Было проведено 4 серии расчетов для смещений 50, 60, 70 и 80 м. При этом для каждого из этих четырех случаев вводились угол падения плоскости разлома 65 град., процентное соотношение «песок / пылеватая глина» 50/50 и поэтапно увеличивалось снижение уровня воды в водоносном горизонте – на 10, 20, 30, 40, 50 и 60 м. На рисунке 10, B показаны результаты расчетов – кривые B1, B2, B3, B4, соответствующие вышеуказанным амплитудам смещений по разлому (в том же порядке). Из рисунка видно, что чем больше смещение по разлому, тем ниже идет кривая, то есть тем больше оседание поверхности земли.

При третьем сценарии при прочих равных условиях оценивалось влияние угла падения плоскости разлома. Было выполнено 3 серии расчетов для углов падения 65, 80 и 50 град. При этом для каждого из этих трех случаев вводились амплитуда смещения по разлому 50 м, процентное соотношение «песок / пылеватая глина» 50/50 и поэтапно увеличивалось снижение уровня воды в водоносном горизонте – на 10, 20, 30, 40, 50 и 60 м. На рисунке 10, C показаны итоги расчетов – кривые С1, С2, С3, соответствующие вышеуказанным углам падения плоскости разлома (в том же порядке). Из рисунка видно, что чем больше угол падения плоскости разлома, тем ниже идет кривая, то есть тем больше осадки поверхности земли.

Кроме того, можно заметить, что амплитуда смещения по разлому влияет на вертикальные смещения меньше, чем угол падения плоскости разлома (см. рис. 10, B, C).

Рис. 10. Результаты параметрического моделирования для территории Анаргири [1]. Влияние на осадки поверхности земли со стороны: А – соотношения мощностей сжимаемых слоев пердикасской (Perdikas) и проастиосской (Proastio) свит, то есть соотношения «песок / пылеватая глина»; B – амплитуды смещения по разлому; C – угла падения плоскости разлома. По оси абсцисс отложены величины снижения уровня грунтовых вод (м), а по оси ординат – оседание поверхности земли (м) в точке В относительно точки А (см. рис. 8) со знаком «минус»

При этом авторы работы [1] подчеркивают, что при всех сценариях параметрического моделирования для территории населенного пункта Анаргири осадки поверхности земли увеличиваются со снижением уровня воды в водоносном горизонте независимо от геотехнических характеристик сжимаемых слоев, амплитуды смещения по разлому и угла падения плоскости разлома (см. рис. 10, A, B, C).

Моделирование для Валтонеры

Участок населенного пункта Валтонера, который находится на расстоянии около 1,5 км от периметра угольного карьера, пересечен вероятным разломом с амплитудой смещения 70 м (см. рис. 2, рис. 11). Этот участок сложен (сверху вниз) аллювиальными отложениями (alluvial deposits), пердикасской (Perdikas) и проастиосской (Proastio) свитами. Однако мощность отложений пердикасской (Perdikas) свиты на юго-восточной стороне (справа от разлома на рисунке 11) существенно больше. Параметры грунтов представлены в таблице 3.

Рис. 11. Упрощенная имитационная 3D-модель грунтовых условий для рассмотренного участка в населенном пункте Валтонера [1]. Мощности слоев указаны в метрах. Геотехнические параметры слоев приведены в таблице 3

Таблица 3. Параметры аллювиальных отложений (alluvial deposits), грунтовых слоев пердикасской (Perdikas) и проастиосской (Proastio) свит на изучаемом участке населенного пункта Валтонера [1]

Результаты численного моделирования вертикальных смещений в грунтовом массиве из-за снижения уровня грунтовых вод в Валтонере представлены на рисунке 12. Эти деформации увеличиваются пропорционально снижению УГВ. Максимальные их значения были отмечены на юго-восточной стороне (справа от разлома), где увеличивается мощность отложений пердикасской (Perdikas) свиты, на которые больше всего влияет падение УГВ.

Рис. 12. Результаты моделирования вертикальных смещений в грунтовом массиве из-за снижения уровня грунтовых вод на рассмотренном участке населенного пункта Валтонера [1]. Цветовая шкала показывает величины вертикальных смещений (м). В таблице в углу рисунка в левом столбце приведены поэтапные величины снижения УГВ (м), а в правом столбце – соответствующие вертикальные смещения точки B относительно точки A (м), отмеченных на модели черным цветом

При понижении уровня грунтовых вод на 10, 20, 30, 40, 50 и 60 м сравнивались величины оседания поверхности земли в точке B относительно точки А (см. рис. 12). Из таблицы в левом углу рисунка 12 видно, что осадки увеличиваются со снижением уровня воды в водоносном горизонте. Они достигают 0,139 м при снижении УГВ на 60 м. Данное значение невелико, а связано это с тем, что отложения пердикасской (Perdikas) свиты, на которые влияет падение УГВ, имеются по обе стороны от разлома.

Для подтверждения достоверности результатов моделирования авторы статьи [1] провели полевые измерения вертикальных смещений в местах разрывов поверхности в Валтонере. Измеренные величины, как и в случае с Анаргири, оказались близкими к полученным при численном анализе методом конечных элементов.

Для рассмотренного участка Валтонеры авторы работы [1] также выполнили параметрическое моделирование по трем сценариям, которые были в целом аналогичны таковым для селения Анаргири.

При первом сценарии при прочих равных условиях оценивалось влияние на осадки поверхности земли со стороны соотношения мощностей сжимаемых отложений пердикасской (Perdikas) и проастиосской (Proastio) свит. Было выполнено 5 серий расчетов для следующих процентных соотношений «песок / пылеватая глина»: 100/0; 70/30; 50/50; 30/70; 0/100. При этом для каждого из этих пяти случаев вводилась амплитуда смещения по разлому 10 м, угол падения плоскости разлома 65 град. и поэтапно увеличивалось снижение уровня воды в водоносном горизонте – на 10, 20, 30, 40, 50 и 60 м. На рисунке 13, А показаны итоги расчетов – кривые А1, А2, А3, А4, А5, соответствующие вышеуказанным соотношениям «песок / пылеватая глина» (в том же порядке). По оси абсцисс отложены величины снижения уровня грунтовых вод (м), по оси ординат – оседание поверхности земли (м) в точке В относительно точки А (см. рис. 12) со знаком «минус». Из рисунка 13, А видно, что чем больше содержание пылевато-глинистых слоев, тем ниже находится график, то есть тем больше оседает поверхность земли.

При втором сценарии при прочих равных условиях изучалось влияние амплитуды смещения по разлому. Было выполнено 7 серий расчетов для смещений 10, 20, 30, 40, 50, 60 и 70 м. При этом для каждого из этих семи случаев вводились угол падения плоскости разлома 65 град., процентное соотношение «песок / пылеватая глина» 50/50 и поэтапно увеличивалось снижение уровня воды в водоносном горизонте – на 10, 20, 30, 40, 50 и 60 м. На рисунке 13, B показаны результаты расчетов – кривые B1, B2, B3, B4, B5, B6, B7, соответствующие вышеуказанным амплитудам смещений (в том же порядке). Из рисунка видно, что чем больше смещение по разлому, тем ниже идет кривая, то есть тем больше оседание поверхности земли.

При третьем сценарии при прочих равных условиях исследовалось влияние угла падения плоскости разлома. Было выполнено 3 серии расчетов для углов падения 65, 80 и 50 град. При этом для каждого из этих трех случаев вводились амплитуда смещения по разлому 20 м, процентное соотношение «песок / пылеватая глина» 50/50 и поэтапно увеличивалось снижение уровня воды в водоносном горизонте – на 10, 20, 30, 40, 50 и 60 м. На рисунке 13, C показаны результаты – кривые С1, С2, С3, соответствующие вышеуказанным углам (в том же порядке). Из рисунка видно, что угол падения плоскости разлома в рассмотренном случае существенно не влияет на вертикальные смещения.

Рис. 13. Результаты параметрического моделирования для территории Валтонеры [1]. Влияние на осадки поверхности земли со стороны: А – соотношения мощностей сжимаемых слоев пердикасской (Perdikas) и проастиосской (Proastio) свит, то есть соотношения «песок / пылеватая глина»; B – амплитуды смещения по разлому; C – угла падения плоскости разлома. По оси абсцисс отложены величины снижения уровня грунтовых вод (м), по оси ординат – оседание поверхности земли (м) в точке В относительно точки А (см. рис. 12) со знаком «минус»

Примечательно, что при всех сценариях моделирования для территории населенного пункта Валтонера (как и для Анаргири) вертикальные деформации увеличиваются пропорционально снижению уровня воды в водоносном горизонте (см. рис. 13, A, B, C).

Моделирование для Фаноса

Участок территории населенного пункта Фанос, который расположен на расстоянии более 4 км от периметра угольного карьера, пересекается разломом F_pxa (см. рис. 2), имеющим амплитуду смещения 80 м. Как видно из рисунка 14, этот участок с поверхности сложен отложениями осыпей, осыпных и аллювиальных конусов выноса (talus cone, scree and alluvial fans), а пердикасская (Perdikas) и проастиосская (Proastio) свиты находятся между этими поверхностными отложениями и неогеновой геологической формацией (Neogene formation) лишь на юго-восточной стороне (справа от разлома), но отсутствуют на северо-западной стороне (слева от разлома). Геотехнические параметры слоев грунта представлены в таблице 4.

Рис. 14. Упрощенная имитационная 3D-модель грунтовых условий для рассмотренного участка в селении Фанос [1]. Мощности слоев указаны в метрах. Геотехнические параметры слоев приведены в таблице 4

Таблица 4. Параметры отложений осыпей, осыпных и аллювиальных конусов выноса (talus cone, scree and alluvial fans), грунтовых слоев пердикасской (Perdikas) и проастиосской (Proastio) свит на изучаемом участке населенного пункта Фанос [1]

Результаты численного моделирования вертикальных смещений в грунтовом массиве из-за снижения уровня грунтовых вод в Фаносе показаны на рисунке 15. Эти деформации увеличиваются пропорционально снижению уровня грунтовых вод. Максимальные смещения наблюдаются на юго-восточной стороне (справа от разлома), где имеются отложения пердикасской (Perdikas) свиты, на которые больше всего влияет снижение УГВ.

Рис. 15. Результаты моделирования вертикальных смещений в грунтовом массиве из-за снижения уровня грунтовых вод на рассмотренном участке населенного пункта Фанос [1]. Цветовая шкала показывает величины вертикальных смещений (м). В таблице в углу рисунка в левом столбце приведены поэтапные величины снижения УГВ (м), а в правом столбце – соответствующие вертикальные смещения точки B относительно точки A (м), отмеченных на модели черным цветом

При понижении уровня воды в водоносном горизонте на 10, 20, 30, 40, 50 и 60 м авторы работы [1] сравнивали величины оседания поверхности земли в точке B относительно точки А (см. рис. 15). Из таблицы в левом углу рисунка 15 видно, что оседание увеличивается со снижением уровня грунтовых вод. Максимальная осадка составила 0,527 м (при снижении УГВ на 60 м). Это значение велико, и это связано с большой разницей между жесткостью неогеновых и четвертичных отложений, а также с тем, что отложения пердикасской (Perdikas) свиты, мощность которых зависит от снижения УГВ, имеются только на юго-восточной стороне (справа от разлома).

Для подтверждения достоверности результатов моделирования Тзапоглу и Лупасакис [1] выполнили полевые измерения вертикальных смещений в местах разрывов поверхности в Фаносе. Измеренные величины, как и для Анаргири и Валтонеры, оказались близкими к полученным при анализе в программе PLAXIS.

Для рассмотренного участка Фаноса авторы статьи [1] тоже провели параметрическое моделирование по трем сценариям, которые были в целом аналогичны таковым для Анаргири и Валтонеры.

При первом сценарии при прочих равных условиях исследовалось влияние на оседание поверхности земли со стороны соотношения мощностей сжимаемых отложений пердикасской (Perdikas) и проастиосской (Proastio) свит. Было проведено 5 серий расчетов для следующих процентных соотношений «песок / пылеватая глина»: 100/0; 70/30; 50/50; 30/70; 0/100. При этом для каждого из этих пяти случаев вводилась амплитуда смещения по разлому 80 м, угол падения плоскости разлома 65 град. и поэтапно увеличивалось снижение уровня воды в водоносном горизонте – на 10, 20, 30, 40, 50 и 60 м. На рисунке 16, А показаны результаты – кривые А1, А2, А3, А4, А5, соответствующие вышеуказанным соотношениям «песок / пылеватая глина» (в том же порядке). По оси абсцисс отложены величины снижения уровня грунтовых вод (м), по оси ординат – оседание поверхности земли (м) в точке В относительно точки А (см. рис. 15) со знаком «минус». Из рисунка 16, А видно, что в целом чем больше содержание пылевато-глинистых слоев, тем ниже находится график, то есть тем больше оседает поверхность земли, но это очевидно только для трех кривых из пяти, в отличие от случаев Анаргири и Валтонеры.

При втором сценарии при прочих равных условиях оценивалось влияние амплитуды смещения по разлому. Было выполнено 3 серии расчетов для смещений 30, 40 и 50 м (так указано у авторов работы [1] на рисунке 16, B, хотя в таблице, которая здесь не приведена, указаны 10, 20 и 30 м. – Ред.). При этом для каждого из этих трех случаев вводились угол падения плоскости разлома 65 град., процентное соотношение «песок / пылеватая глина» 50/50 и поэтапно увеличивалось снижение УГВ – на 10, 20, 30, 40, 50 и 60 м. На рисунке 16, B показаны результаты расчетов – кривые B1, B2, B3, соответствующие вышеуказанным амплитудам смещений (в том же порядке). Из рисунка видно, что чем больше смещение по разлому, тем ниже идет кривая, то есть тем больше оседание поверхности земли.

При третьем сценарии при прочих равных условиях рассматривалось влияние угла падения плоскости разлома. Было выполнено 3 серии расчетов для углов падения 65, 80 и 50 град. При этом для каждого из этих трех случаев вводились амплитуда смещения по разлому 80 м, процентное соотношение «песок / пылеватая глина» 50/50 и поэтапно увеличивалось снижение УГВ – на 10, 20, 30, 40, 50 и 60 м. На рисунке 16, C показаны результаты – кривые С1, С2, С3, соответствующие вышеуказанным углам (в том же порядке). Из рисунка видно, что угол падения плоскости разлома в рассмотренном случае не влияет на вертикальные смещения поверхности земли.

Рис. 16. Результаты параметрического моделирования для территории Фаноса [1]. Влияние на осадки поверхности земли со стороны: А – соотношения мощностей сжимаемых слоев пердикасской (Perdikas) и проастиосской (Proastio) свит, то есть соотношения «песок / пылеватая глина»; B – амплитуды смещения по разлому; C – угла падения плоскости разлома. По оси абсцисс отложены величины снижения уровня грунтовых вод (м), по оси ординат – оседание поверхности земли (м) в точке В относительно точки А (см. рис. 15) со знаком «минус»

При этом авторы работы [1] особо отмечают, что при всех сценариях моделирования для территории населенного пункта Фанос (как и для Анаргири и Валтонеры) вертикальные деформации увеличиваются пропорционально снижению уровня воды в водоносном горизонте (см. рис. 16, A, B, C).

ЗАКЛЮЧЕНИЕ

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

Изучение авторами работы [1] четвертичных отложений на территории угленосного суббассейна Аминдео по литературным данным и результатам полевых и лабораторных измерений показало, что пылевато-глинистые грунты с высокими значениями индекса компрессии С_с, что особенно касается пердикасской (Perdikas) свиты, подвержены уплотнению в зонах, где происходит снижение уровня воды в водоносном горизонте.

Исследования показали, что распределение по изучаемой территории тектонических линий (см. рис. 2) согласуется с разрывами поверхности, имеющими вертикальные смещения. Это говорит о том, что на оседание поверхности земли помимо механических свойств геоматериалов влияют также тектонические разломы.

Оценка процесса изменений УГВ привела авторов статьи [1] к выводу, что дренаж путем откачки воды через скважины при разработке угольного карьера Аминдео, а также забор воды из скважин системы орошения для полива сельскохозяйственных культур в течение последних десятилетий привели к весьма чрезмерной эксплуатации водоносного горизонта. Снижение уровня грунтовых вод распространяется на расстояние порядка 4 км от карьера и достигает максимальных значений 70–100 м рядом с ним, образуя широкую депрессионную воронку, в пределах которой оказался ряд населенных пунктов, в том числе Анаргири, Валтонера и Фанос. Более того, в течение последних десятилетий УГВ имеет тенденцию к постоянному снижению (из-за непрерывного потока грунтовых вод в сторону карьера) с сезонными колебаниями (видимо, связанными с работой скважин систем орошения сельскохозяйственных земель).

Анализ методом конечных элементов в программном комплексе PLAXIS 3D (версии 2016 года) показал, что степень влияния основных причинных факторов на величины осадок поверхности земли на разных участках рассматриваемой территории различна.

Численное моделирование для населенного пункта Анаргири показало большие вертикальные деформации из-за того, что его территория пересекается разломом со смещением около 50 м, нарушающим геологическое строение. Пердикасская (Perdikas) свита отсутствует с юго-восточной стороны рассмотренного участка относительно разлома. Параметрическое моделирование по трем сценариям показало, что на осадки сильно влияет соотношение мощностей сжимаемых отложений пердикасской (Perdikas) и проастиосской (Proastio) свит, а точнее доля (по мощности) слоев пылеватых глин в разрезе. Было также выявлено, что угол падения плоскости разлома влияет на вертикальные смещения больше, чем амплитуда смещения по разлому.

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

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

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

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

• снижение уровня грунтовых вод;
• расположение и сжимаемость грунтовых слоев;
• тектоническое строение и амплитуды смещений по разломам.

Тзапоглу и Лупасакис [1] с удовлетворением отмечают, что сочетание классических инженерно-геологических лабораторных и полевых исследований с современными методиками моделирования на основе метода конечных элементов помогло им не только понять механизм оседания поверхности земли, но и определить степени влияния на него основных факторов. При этом авторы работы [1] подчеркивают, что отдельные анализы с помощью МКЭ в подобных случаях следует проводить прежде всего для участков территорий, где изменяются гидрогеологические и тектонические условия. Полученная таким способом информация будет полезна для органов и лиц, принимающих решения, дав им возможность вовремя принимать необходимые меры для безопасного развития территорий.

Источники

1. Tzampoglou P., Loupasakis C. Numerical simulation of the factors causing land subsidence due to overexploitation of the aquifer in the Amyntaio open coal mine, Greece // HydroResearch. Vol. 1. Elsevier B.V.: June 2019. P. 8–24. URL: sciencedirect.com/science/article/pii/S2589757819300034.

Список литературы, использованной авторами статьи [1]

Abidin, H.Z., Andreas, H., Gumilar, I.,Wangsaatmaja, S., Fukuda, Y., Deguchi, T., 2009. Land subsidence and groundwater extraction in Bandung Basin, Indonesia. IAHS Publ. 20, 145.

Abidin, H.Z., Andreas, H., Gumilar, I., Yuwono, B.D., Murdohardono, D., Supriyadi, S., 2015. On Integration of Geodetic Observation Results for Assessment of Land Subsidence Hazard Risk in Urban Areas of Indonesia. IAG 150 Years. vol. 143 pp. 435–442.

Al Heib, M., Duval, C., Theoleyre, F., Watelet, J.-M., Gombert, P., 2015. Analysis of the historical collapse of an abandoned underground chalk mine in 1961 in Clamart (Paris, France). Bull. Eng. Geol. Environ. 74 (3), 1001–1018.

Antoniadis, P., Mavridou, E., Gentzis, T., 2005. The Notio Pedio (Southern Field) lignite deposit in the Ptolemaida Basin, Greece: depositional conditions as revealed through petrography. Energy Sources 27 (12), 1117–1131.

Bell, J.W., Amelung, F., Ramelli, A.R., Blewitt, G., 2002. Land subsidence in Las Vegas, Nevada, 1935–2000: newgeodetic data showevolution, revised spatial patterns, and reduced rates. Environ. Eng. Geosci. 8 (3), 155–174.

Choi, J.-K., Kim, K.-D., Lee, S.,Won, J.-S., 2010. Application of a fuzzy operator to susceptibility estimations of coal mine subsidence in Taebaek City, Korea. Environ. Earth Sci. 59 (5), 1009–1022.

Conway, B.D., 2016. Land subsidence and earth fissures in south-central and southern Arizona, USA. Hydrogeol. J. 24 (3), 649–655.

Deverel, S.J., Rojstaczer, S., 1996. Subsidence of agricultural lands in the Sacramento-San Joaquin Delta, California: role of aqueous and gaseous carbon fluxes. Water Resour. Res. 32 (8), 2359–2367.

Dhmaras, K., Georgiadhs, M., 2002. Lifting and Waterproofing Study of the Existing embankment of Lake Cheimaditida Florina, Geognosi A.E., Thessaloniki.

Dimitrakopoulos, D., 2001. Hydrogeological conditioning of Amyndeon mine. Problems during exploitation and overcoming them. Phd Thesis. National Technical University of Athens, Athens.

Doukissa, K., 2010. Geotechnical Study of the Bridge in the Central Trench of Aanargiroi Area, Municipality of Aetos, Thessaloniki.

Fakhri, F., Kalliola, R., 2015. Monitoring ground deformation in the settlement of Larissa in Central Greece by implementing SAR interferometry. Nat. Hazards 78 (2), 1429–1445.

Galloway, D.L., Burbey, T.J., 2011. Review: regional land subsidence accompanying groundwater extraction. Hydrogeol. J. 19 (8), 1459–1486.

Galloway, D.L., Jones, D.R., Ingebritsen, S.E., 1999. Land subsidence in the United States. US Geological Survey. Reston, Virginia.

Ganas, A., Salvi, S., Atzori, S., Tolomei, C., 2006. Ground deformation in Thessaly, Central Greece, retrieved from differential interferometric analysis of ERS-SAR data. 11th International Symposiumon Natural and Human Induced Hazards & 2ndWorkshop on Earthquake Prediction, Patras, Greece. 41, pp. 22–25.

Ganas, A., Oikonomou, I., Tsimi, A., 2015. NOAfaults: a digital database for active faults in Greece. Bull. Geol. Soc. Greece 47, 518–530.

Geoerevnitiki, E., 1980. AHS Amyntaiou, Preliminary Geotechnical Investigation. PPC Athens, Athens.

Gong, S., LI, P., Wu, J., Li, X., 2009. Analysis of groundwater seepage field and stratum stress field of land subsidence in Shanghai. J. Liaoning Technical Univ. 29, 232–235 Natural Science).

Gongyu, L., Wanfang, Z., 1999. Sinkholes in karst mining areas in China and some methods of prevention. Eng. Geol. 52 (1), 45–50.

Guerrero, J., Gutiйrrez, F., Bonachea, J., Lucha, P., 2008. A sinkhole susceptibility zonation based on paleokarst analysis along a stretch of theMadrid–Barcelona high-speed railway built over gypsum-and salt-bearing evaporites (NE Spain). Eng. Geol. 102 (1), 62–73.

Herrera, G., Tomбs, R., Monells, D., Centolanza, G., Mallorquн, J., Vicente, F., Navarro, V., Sanchez, J., Sanabria, M., Cano, M., Mulas, J., 2010. Analysis of subsidence using TerraSAR-X data: Murcia case study. Eng. Geol. 116, 284–295.

Huang, B., Shu, L., Yang, Y., 2012. Groundwater overexploitation causing land subsidence: hazard risk assessment using field observation and spatial modelling. Water Resour. Manag. 26 (14), 4225–4239.

I.G.M.E, 1997. Geological Map of Greece, Scale 1:50.000, Ptolemaida Sheet. Institute of Geology and Mineral Exploration, Athens.

Kaitantzian, A., Loupasakis, C., Rozos, D., 2014. Assessment of geo-hazards triggered by both natural events and human activities in rapidly urbanized areas. In: Lollino, G., et al. (Eds.), Engineering Geology for Society and Territory. Torino, Italy, pp. 675–679.

Kontogianni, V., Pytharouli, S., Stiros, S., 2007. Ground subsidence, Quaternary faults and vulnerability of utilities and transportation networks in Thessaly, Greece. Environ. Geol. 52 (6), 1085–1095.

Koukouzas, C., Kotis, T., Ploumidis, M., Metaxas, A., 1979. Coal exploration of Anargiri- Amynteon area. Mineral Deposit Research. Institute of Geology and Mineral Exploration, Athens.

Koukouzas, C., Kotis, T., Ploumidis,M.,Metaxas, A., Dimitriou, D., 1983. Coal Exploration of Komninon-Ptolemaidas Area. Mineral Deposit Research. Institute of Geology and Mineral Exploration, Athens.

Li, Q.,Wang, H., 2006. A study on land subsidence in Shanghai. Geol. J. China Univ. 12 (2), 169–178.

Li, Y., Gong, H., Zhu, L., Li, X., 2017. Measuring spatiotemporal features of land subsidence, groundwater drawdown, and compressible layer thickness in Beijing Plain, China. Water. 9 (1), 64.

Loupasakis, C., 2006. Study of the geotechnical conditions of the Amintaio coalmine slopes close to the Anargiri village, Aetos Municipality, Florina Prefecture, Greece. Unpublished Report. Institute of Geology and Mineral Exploration, Athens.

Loupasakis, C., 2010. Geotechnical works in the Anargirпi village, Aetos Municipality, Florina Prefecture, Greece. Unpublished Report. Institute of Geology and Mineral Exploration, Athens.

Loupasakis, C., Rozos, D., 2009. Finite-element simulation of land subsidence induced by water pumping in Kalochori village, Greece. Q. J. Eng. Geol. Hydrogeol. 42 (3), 369–382.

Loupasakis, C., Angelitsa, V., Rozos, D., Spanou, N., 2014. Mining geohazards—land subsidence caused by the dewatering of opencast coal mines: the case study of the Amyntaio coal mine, Florina, Greece. Nat. Hazards 70 (1), 675–691.

Mavridou, E., Antoniadis, P., Khanaqa, P., Riegel,W., Gentzis, T., 2003. Paleoenvironmental interpretation of the Amynteon–Ptolemaida lignite deposit in northern Greece based on its petrographic composition. Int. J. Coal Geol. 56 (3), 253–268.

Metaxas, A., Karageorgiou, D., Varvarousis, G., Kotis, T., Ploumidis, M., Papanikolaou, G., 2007. Geological evolution-stratigraphy of Florina, Ptolemaida, Kozani and Saradaporo graben. Bull. Geol. Soc. Greece 40 (1), 161–172.

Motagh, M., Walter, T.R., Sharifi, M.A., Fielding, E., Schenk, A., Anderssohn, J., Zschau, J., 2008. Land subsidence in Iran caused by widespread water reservoir overexploitation. Geophys. Res. Lett. 35 (16).

Mountrakis, D., 1985. Geology of Greece, University studio Press, Thessaloniki.

Mountrakis, D., Pavlides, S., Zouros, N., Astaras, T., Chatzipetros, A., 1998. Seismic fault geometry and kinematics of the 13May 1995Western Macedonia (Greece) earthquake. J. Geodyn. 26 (2), 175–196.

Mountrakis, D., Tranos, M., Papazachos, C., Thomaidou, E., Karagianni, E., Vamvakaris, D., 2006. Neotectonic and seismological data concerning major active faults, and the stress regimes of Northern Greece. Geol. Soc. Lond., Spec. Publ. 260 (1), 649–670.

Nieuwenhuis, H., Schokking, F., 1997. Land subsidence in drained peat areas of the Province of Friesland, The Netherlands. Q. J. Eng. Geol. Hydrogeol. 30 (1), 37–48.

Pacheco-Martнnez, J., Hernandez-Marнn, M., Burbey, T.J., Gonzбlez-Cervantes, N., Ortнz- Lozano, J.Б., Zermeсo-De-Leon, M.E., Solнs-Pinto, A., 2013. Land subsidence and ground failure associated to groundwater exploitation in the Aguascalientes Valley, Mexico. Eng. Geol. 164, 172–186.

Parcharidis, I., Foumelis, M., Katsafados, P., 2011. Seasonal Ground Deformation Monitoring Over Southern Larissa Plain (Central Greece) by SAR Interferometry. Advances in the Research of Aquatic Environment. vol. 2 pp. 497–504.

Pavlides, S., 1985. Neotectonic evolution of the Florina-Vegoritis-Ptolemais basin (W. Macedonia, Greece). Phd Thesis. Aristotle University of Thessaloniki, Thessaloniki.

Pavlides, S., Mountrakis, D., 1985. Neotectonics of the Florina-Vegoritis-Ptolemais Neogene Basin (NWGreece): An Example of Extensional Tectonics of the Greater Aegean Area. Laboratoire de gйologie de l'Universitй. pp. 311–327.

Pavlides, S., Mountrakis, D., 1987. Extensional tectonics of northwestern Macedonia, Greece, since the late Miocene. J. Struct. Geol. 9 (4), 385–392.

Pavlides, S., Simeakis, K., 1988. Neotectonics and active tectonics in lowseismicity areas of Greece: Vegoritis (NW Macedonia) and Melos isl. Complex–Comparison. Ann. Gйol. Pays Hellйn. 33 (191–176).

Raspini, F., Loupasakis, C., Rozos, D., Adam, N., Moretti, S., 2014. Ground subsidence phenomena in the Delta municipality region (Northern Greece): geotechnical modeling and validation with persistent Scatterer interferometry. Int. J. Appl. Earth Obs. Geoinf. 28, 78–89.

Raspini, F., Bianchini, S., Moretti, S., Loupasakis, C., Rozos, D., Duro, J., Garcia, M., 2016. Advanced interpretation of interferometric SAR data to detect, monitor and model ground subsidence: outcomes from the ESA-GMES Terrafirma project. Nat. Hazards 83 (1), 155–181.

Rothenburg, L., Obah, A., El Baruni, S., 1995. Horizontal ground movements due to water abstraction and formation of earth fissures. International Association of Hydrological Sciences. IAHS Publ. 234, 239–249.

Rozos, D., Sideri, D., Loupasakis, C., Apostolidhs, E., 2010. Land subsidence due to excessive ground water withdrawal. A case study from Stavros-Farsala site, west Thessaly Greece. Proceedings of the 12th International Congress of the Geological Society of Greece. Bull. Geol. Soc. Greece 43 (4), 1850–1857.

Soulios, G., Tsapanos, T., Voudouris, K., Kaklis, T., Mattas, C., Sotiriadis, M., 2011. Ruptures on surface and buildings due to land subsidence in Anargyri village (Florina Prefecture, Macedonia). Environ. Earth Sci. 5, 505–512.

Spyropoulos, N., 1992. The geological structure of Pelagonian zone in the Mount of Askio D. Macedonia. Phd Thesis. Aristotle University of Thessaloniki, Thessaloniki.

Stamos, A., Giannoulopoulos, P., 2010. Hydrogeological report, Geotechnical work to the Anargiri village, Aetos Municipality, Florina Prefecture, Greece. Institute of Geology and Mineral Exploration, Athens.

Strzalkowski, P., Tomiczek, K., 2015. Analytical and numerical method assessing the risk of sinkholes formation in mining areas. Int. J. Min. Sci. Technol. 25 (1), 85–89.

Svigkas, N., Papoutsis, I., Loupasakis, C., Tsangaratos, P., Kiratzi, A., Kontoes, C., 2016. Land subsidence rebound detected via multi-temporal InSAR and ground truth data in Kalochori and Sindos regions, Northern Greece. Eng. Geol. 209, 175–186.

Tosi, L., Lio, C., Teatini, P., Strozzi, T., 2018. Land subsidence in coastal environments: knowledge advance in the Venice coastland by TerraSAR-X PSI. Remote Sens. 10 (8), 1191.

Tsourlos, P., 2015. Research for the appearance of surface cracks in the Valtonera village, Amyntaio Мunicipality, Florina Region. School of Geology of the Aristotle University of Thessaloniki Aristotelian University of Thessaloniki, Thessaloniki.

Tsourlos, P., Papazaxos, C., Bargmeziw, G., Koutalou, B., 2007. Investigation and suggestions for the problem created by the fault which as appeared in municipalities of Fanos and Ksino Nero. School of Geology of the Aristotelian University of Thessaloniki, Thessaloniki.

Tzampoglou, P., Loupasakis, C., 2016. New data regarding the ground water level changes at the Amyntaio basin-Florina Prefecture, Greece. Bull. Geol. Soc. Greece 50 (2), 1006–1015.

Tzampoglou, P., Loupasakis, C., 2017a. Evaluating geological and geotechnical data for the study of land subsidence phenomena at the perimeter of the Amyntaio coalmine, Greece. Int. J. Min. Sci. Technol. 28 (4), 61–612.

Tzampoglou, P., Loupasakis, C., 2017b. Land Subsidence Susceptibility and Hazard Mapping: The Case of Amyntaio Basin, Greece, International Society for Optics and Photonics. p. 104441L.

Tzampoglou, P., Loupasakis, C., 2017c. Mining geohazards susceptibility and risk mapping: the case of the Amyntaio open-pit coalmine,WestMacedonia, Greece. Environ. Earth Sci. 76 (15), 542.

Tzampoglou, P., Loupasakis, C., 2017d. Updated Ground Water Piezometry Data of the Amyntaio Sub-Basin and their Effect to the Manifestation of the Land Subsidence Phenomena. 11th International Hydrogeological Congress of Greece, Athens.

Vetoulis, D., 1951. Lignite Basin of Amyntaio- Vegora, Ed. Geological IGEF Recognition.

Wanfang, Z., 1997. The formation of sinkholes in karst mining areas in China and some methods of prevention. Environ. Geol. 31 (1–2), 50–58.

Wang, H.G., 2006. Current situation and tendency of land subsidence in Tianjin. Master Dissertation. China University Geoscience.

Wei, J.S., Xu, J., Lu, Y., Yi, C.R., 2012. Analysis on land subsidence funnel in Wangqingtuo area of Tianjin. Ground Water 34 (5), 49–51.

Yiping,W., 2004. Study on land subsidence caused by overexploiting groundwater in Beijing. Site Investig. Sci. Technol. 5, 46–49.

Zhang, X., 2000. On establishment of Suzhou earth subsidence gray forecast model. J. Suzhou Institute of Urban Construction and Environmental Protection. 13 (4), 53–57.

Zhang, Z.Y., 2007. Analysis on present situation and development trend of groundwater descent funnel in Cangzhou City. Ground Water 29 (1), 50–53.