Introduction

최근 국내에서 발생하는 지반침하(이른바 싱크홀, sinkhole)의 주요 원인은 자연적 요인보다 인위적 요인에 기인하는 것으로 보고되고 있다. 국토안전관리원(Korea Authority of Land & Infrastructure Safety, KALIS) 및 보도 자료에 따르면, 최근 10년간 지반침하 사례의 약 절반이 상·하수도관 손상에 따른 누수에서 비롯된 것으로 집계된다. 예컨대 2014 - 2023년 통계에서는 하수관 손상 42.01%, 상수관 손상 12.61%로 합계 54.62%가 보고되었고(Kim, 2025), 2019 - 2023년 자료에서도 전체 957건 중 하수관 446건, 상수관 39건으로 합계가 약 50%를 상회한다(KALIS, 2024; Water Journal, 2024). 장기간의 부식·균열·비틀림으로 인한 누수는 관로 주변 토사를 유실시켜 공동(void)을 형성하고, 결과적으로 지반침하로 이어진다. 특히 국내 상·하수도관 중 약 40%가 설치 후 20년 이상 경과한 노후관으로 분류되어 향후 위험 증대를 시사한다(KALIS, 2024).

해외 사례도 동일한 맥락을 시사한다. 프랑스 오를레앙에서는 2016년 홍수 기간 불과 열흘 사이 12건의 싱크홀이 집중 발생해, 집중호우와 지하수위 변동이 공동 발달과 붕괴를 가속화할 수 있음을 보여주었다(Luu et al., 2019). 일본 후쿠오카 JR 하카다역 전면부에서는 2016년 지하철 연장 공사 중 지하수 유입으로 도로가 대규모로 함몰되었고, 이후 레이더 기반 공동 조사의 정례화와 하수관 갱신이 강화되었다(Ichii, 2025). 2025년 1월 사이타마현 야시오시의 직경 약 5 m, 깊이 약 10 m 도로 함몰 사례에서는 조사위원회가 노후 하수관 부식을 주요 원인으로 지목했으며(NHK, 2025), 국토교통성은 사회적 영향이 큰 구간의 관로 갱신을 특별 보조 대상으로 지정하고 노후 하수관로에 대한 전국적 집중 조사를 추진하고 있다.

이와 관련된 지반침하 연구는 주로 현장조사와 실내실험에 의존해 왔으며, 집중호우·지하수위 변동·지하시설물 열화와 같은 외적 교란에 따른 공동 발생과 붕괴 메커니즘을 사건별로 기록하는 접근이 다수를 차지한다(Luu et al., 2019). 제한된 현장자료 내에서 다양한 물성치·수리조건을 체계적으로 교차 검증하려면, 물리실험만으로는 시간·비용·안전에 대한 한계가 크므로 수치모형을 병행하는 것이 합리적이다(Zhang et al., 2016). 특히 지반침하와 같이 국부적 대변형, 공극수 유동에 의한 강한 비선형성, 접촉·분리·균열이 수반되는 문제에는 재료점법(material point method, MPM)의 적용이 적절하다. MPM은 입자(재료점)-오일러 격자 하이브리드 기법으로, 라그랑주법의 물성 이력 추적 용이성과 오일러법의 격자 왜곡 회피 장점을 결합해 대변형과 파괴, 다상 상호작용을 견고하게 계산할 수 있도록 정식화되어 있다(Zhang et al., 2016).

이상의 국내외 현황을 살펴보면 지하수위 변동, 집중호우, 관로 노후화·열화가 복합적으로 작용해 도시 환경에서 싱크홀이 유발됨을 뚜렷이 시사한다. 그럼에도 관망 주변에서의 지하 공동 발생을 기점으로 표층 붕괴에 이르는 시간적 변화를 고해상도 수치해석으로 연속적으로 재현한 연구는 제한적이다. 본 연구는 이 공백을 보완하고자, 관망 주변 지하공동 발생 후 붕괴 개시를 가정한 수치모의를 수행해 붕괴의 시간적 변화를 분석하였다. 이를 위해 MPM의 기본 체계(Zhang et al., 2016)를 참조하여 코드를 구현하고, 병렬처리를 적용하였다. 계산도메인의 토층의 경우 최상부 지표면은 도로가 위치한다고 가정하여 도로 포장 재료인 아스팔트 층으로, 그 하부에는 일반적인 사질토 층이 존재한다고 가정하여 재료구성모델과 경계조건을 설정하였다. 본 접근은 공동 주변의 대변형과 비선형 거동을 시간사상으로 해석하여, 향후 도시 지반 안정성 평가와 예방정책 수립을 위한 기계적 근거를 제공한다.

Materials and Methods

MPM model

본 연구는 연속체 기반 입자법인 MPM의 기본 알고리즘을 채택하였다. MPM은 라그랑주 입자와 오일러 배경격자를 결합해 대변형과 파괴, 분절화를 안정적으로 모의하며, 약식(weak form) 기반으로 구성방정식의 이력 의존성을 입자에서 해석하는 특징을 가진다(Zhang et al., 2016). MPM의 이론적 구성(갱신 라그랑주 정식화, 운동방정식 약형, 입자-격자 사상 및 재사상, 접촉 및 비압축성 처리 등)과 구현 지침은 본 연구에서 활용한 MPM 모형의 매뉴얼(Zhang et al., 2016)에 상세히 제시되어 있으므로(코드 구조와 예제 등) 세부 내용은 이를 참고할 수 있다. 본 연구에서는 위 기본 코드에 계산 성능 향상을 위해 OpenMP 병렬처리를 적용하고 Paraview 가시화를 위한 Visualization Toolkit unstructured grid (VTU) 바이너리 파일 출력을 추가하였다.

Simulation conditions and scenarios

Table 1처럼 본 연구의 모의는 3차원 도메인(20 × 15 × 2 m³)에 기반하며 격자간격(정육면체)은 0.5 m이다. 매 시간간격 마다 Courant-Friedrichs-Lewy (CFL) 조건과 dashpot-spring에 따른 가변 dt로 결정된다. 총 해석 시간은 5 s이며 경계조건은 도메인 내부에서만 거동이 발생하고 외부유출이 없는 폐쇄 경계조건이다. 여기서, 각 경계면에 접촉된 입자들은 입자와 동일한 물성치의 입자가 경계면에 작용하고 있는 것으로 가정한다.

본 연구에서는 수치모의에 필요한 입력 파라미터를 국내 설계기준(KDS)에 근거하여 설정하였다(Table 2). 사질토에 대해서는 ｢KDS 11 00 00 지반설계기준(MOLIT, 2021a)｣에서 명시하듯이 전단강도정수(φ′, c)는 현장시험 또는 상관식으로 산정하도록 규정되어 있으며, 고정된 표준값은 제시되어 있지 않다. 이에 따라 국내 일반 설계에서 널리 활용되는 중밀-조밀 사질토의 내부마찰각 범위(약 30 - 35°)를 참고하여, 본 연구에서는 φ′ = 33° (≈0.58 rad)를 채택하였다. 점착력은 배제(c = 0), 팽창각 또한 무시(ψ = 0)하였다. 변형계수는 정착 및 침하 검토 시 활용되는 변형계수(Es)에 근거하여 보수적으로 약 70 MPa를 적용하였다.

아스팔트 혼합물의 경우, ｢KDS 44 50 05 도로포장설계기준(MOLIT, 2021b)｣은 동탄성계수(E*)와 같은 성능 지표를 통해 재료 특성을 반영하도록 규정하고 있어 Mohr-Coulomb 파라미터(c-φ-ψ)는 명시되어 있지 않다. 그러나 본 연구에서 사용한 MPM 해석 구조상 Mohr-Coulomb 파라미터 입력이 요구되므로, 준정적 20℃ 조건에서 기존 문헌(Huang, 2004; KSA, 2014; MOLIT, 2020)에서 제시된 아스팔트 콘크리트의 물성 범위를 근거로 등가 파라미터를 추정하였다. 구체적으로, 탄성계수는 E = 3.0 GPa (일반적으로 2.0 - 6.0 GPa 범위), 포아송비는 ν = 0.35 (0.30 - 0.40 범위), 내부마찰각은 φ = 20° (≈0.35 rad, 보고치 18 - 25° 범위), 유효 점착력은 c = 0.5 MPa (0.3 - 1.0 MPa 범위), 인장강도는 ft = 0.7 MPa (0.5 - 1.0 MPa 범위)로 설정하였다. 이 값은 KDS의 동탄성 계수 개념을 정역학적 파라미터 체계로 환산·근사한 모의용 계수로, 실제 설계용 강도 정수라기보다는 수치해석적 단순화를 위한 입력치라는 점을 명확히 한다.

이와 같이 아스팔트와 사질토는 KDS 절차 준수 전제 하의 보수적 설계치로, 아스팔트는 기준 부재 항목에 대한 모의용 추정치로 제시하였다.

본 연구에서는 공동이 존재한다는 가정을 기반으로 공동의 직경을 변화시켜 지반 거동을 비교하기 위해 네 가지 수치모의 시나리오를 설정하였다(Table 3; Fig. 1). Run1 (직경: 10 m)에서 Run4 (직경: 4 m)로 갈수록 공동 직경이 작아지고 물질점 수는 증가하여, 지반 내 공극 확대에 따른 상부 지반 지지력 상실 과정을 단계적으로 모의하였다.

Fig. 1.

Computational domain.

Results and Discussion

Time changes in simulation result

본 절에서는 공동이 가장 큰 경우인 Run1 (직경 10 m)의 시간의 경과에 따른 지반 거동을 Fig. 2에 나타냈다. 초기 상태(Fig. 2A)에서는 상부 지반이 안정적으로 유지되었으나, 1 s 시점(Fig. 2B)부터 공동 천단부에 응력이 집중되면서 표층 입자가 점진적으로 낙하하기 시작하였다. 이후 2 s (Fig. 2C)에는 공동 상부의 지지력이 상실되며 대규모 입자 붕괴가 발생하였고, 3 s (Fig. 2D)에서는 낙하된 입자가 공동 내부를 채우는 동시에 상부 지반의 추가적인 균열이 확산되었다. 4 s 시점(Fig. 2E)에서는 함몰부의 상단이 확대되면서 불안정한 균형 상태가 형성되었으며, 최종적으로 5 s (Fig. 2F)에는 지표면이 급격히 함몰하여 전형적인 싱크홀 형태가 나타났다. 이 과정은 실제 대규모 도로 함몰 사례에서 보고된 단시간 내 붕괴 및 침하 현상과 유사한 메커니즘을 재현한 것으로 볼 수 있다. 5 s 시점에서는 입자의 위치에 따라 붕괴 양상이 뚜렷하게 구분되었다. 도로 가장자리 경계부의 입자는 점착력의 영향으로 일정 시간 동안 잔류하며 안정적인 모습을 보였으나, 중앙부 상부에 위치한 두 줄의 입자는 즉시 낙하하여 공동 내부로 빠르게 충전되었다. 이러한 결과는 싱크홀 발생 시 중앙부는 지지력이 거의 없어 급격한 함몰이 일어나는 반면, 가장자리는 상대적으로 지반의 결합력이 유지되어 부분적으로 잔류할 수 있음을 보여준다. 따라서 실제 도로 함몰 사고에서도 중앙부 급격 침하와 주변부의 지연 붕괴가 동시에 관찰될 수 있음을 시사한다.

Fig. 2.

Time changes in sinkhole behavior of Run1.

Run1의 수치모의 결과는 공동 직경 10 m 조건에서 중앙부 입자의 즉각적인 낙하와 도로 가장자리 경계부 입자의 점착력에 의한 잔류라는 붕괴 양상을 보여주었다. 이러한 결과는 실제 대규모 도로 함몰 사례와 유사하다. 2016년 후쿠오카 JR 하카다역 앞 사고는 지하철 연장 공사 중 터널 상부 지반이 붕괴하면서 길이 약 30 m, 폭 약 27 m의 싱크홀이 단시간에 발생한 사례이다. 이때 중앙부는 급격히 붕괴하였으나 주변부는 점차 안정화되는 양상이 보고되었는데(Katz, 2016; Wired, 2016), 이는 본 연구에서 관찰된 “중앙부 급속 함몰-경계부 지연 붕괴” 패턴과 일치한다. 2025년 사이타마현 야시오시에서는 노후 하수도관의 부식으로 공동이 형성되어 직경 약 10 m, 깊이 약 10 m의 도로 함몰이 발생하였고, 이후 최대 약 40 m까지 확장되었다(Gayle, 2025a; 2025b; Hashimoto, 2025). 조사 결과 내측 무근 콘크리트는 완전히 손실된 반면 외측 철근부는 일부 잔류하고 있었으며, 중앙부 붕괴가 우선적으로 발생한 점은 본 연구 결과와 구조적으로 대응한다.

따라서 본 연구의 수치모의는 실제 사고 사례에서 확인된 전형적 붕괴 메커니즘을 재현하고 있으며, 향후 싱크홀 위험평가 시 공동 중심부의 지지력 상실뿐만 아니라 경계부의 잔류 및 지연 붕괴 특성까지 고려할 필요가 있다.

Final patterns in simulations

Fig. 3은 공동 직경 조건에 따른 Run2, Run3, Run4의 최종 시점 함몰 형상을 비교한 것이다. Run2 (Fig. 3A)의 경우 함몰부의 깊이는 약 8.42 m로 측정되었으며, 붕괴가 공동 중심부에서 시작하여 주변부까지 확산되면서 지표면이 광범위하게 침하하는 양상을 보였다. Run3 (Fig. 3B)에서는 함몰부 깊이가 약 4.39 m로 Run2보다 얕게 형성되었으며, 붕괴가 중앙부에 집중되는 형태를 나타냈다. 이는 공동 직경이 확대됨에 따라 상부 토층이 부분적으로 아치(arch) 구조를 형성하여 붕괴가 국부적으로 제한된 결과로 해석된다. 반면 Run4 (Fig. 3C)에서는 지표면의 침하 깊이가 약 0.41 m에 불과하였으나, 내부에서는 약 1.86 m의 수직 방향 변형이 발생하였다. 이는 지표면이 점착력 및 경계부 지지력에 의해 단기적으로 유지되었음에도 불구하고, 내부에서는 이미 붕괴가 진전되고 있었음을 보여준다.

Run2에서 Run4로 갈수록 공동 직경은 증가하였음에도 불구하고 지표면의 침하 폭은 오히려 감소하는 경향을 보였다. 즉, 공동 규모가 커질수록 함몰은 지표 전면에서 즉각적으로 발현되기보다는 중앙부 집중형 혹은 내부 잠재 붕괴 형태로 전환됨을 확인할 수 있었다. 이러한 결과는 실제 후쿠오카(2016) 및 사이타마(2025) 도로 함몰 사례와 유사하다. 후쿠오카 사고에서는 공동 확장 이후 지표면이 단기간에 광범위하게 붕괴한 반면(Katz, 2016; Wired, 2016), 사이타마 사고에서는 노후 하수관 부식으로 대규모 공동이 형성되었음에도 불구하고 초기 함몰은 중앙부에서 국부적으로 시작되었으며, 가장자리 구간은 잔류하는 양상을 보였다(Gayle, 2025a; 2025b; Hashimoto, 2025). 이는 본 연구의 수치모의 결과가 실제 대규모 도로 함몰 메커니즘—즉, 중앙부 급속 붕괴와 경계부 지연 붕괴—를 재현하는 데 유효함을 시사한다.

Fig. 3.

Final patterns in Run2 – 4.

Residual voids and secondary failure risk

Run4의 경우 지표면 침하 깊이는 약 0.41 m로 제한적이었으나, 내부에서는 약 1.86 m의 수직 변형이 발생하였다. 특히 주목할 점은 공동이 최종 단계에서도 완전히 충전되지 않은 상태로 남아 있다는 것이다. 이는 향후 지반 거동에 있어 두 가지 불안정 시나리오를 내포한다. 첫째, 공동 잔류 부위가 장기간 유지될 경우 추가적인 토사 유실이나 지하수 유입 등에 의해 공동이 더 확대될 가능성이 존재한다. 둘째, 상부 지반 침하가 발생하여 해당 구간을 단순히 보수할 경우, 보수층의 하중에 의해 잔류 공동이 압밀·붕괴되면서 2차 침하가 재차 발생할 수 있다.

이러한 특성은 실제 사이타마(2025) 사고에서도 관찰된 바 있으며, 당시 내측 무근 콘크리트가 손실된 상태에서 외측 철근부가 일시적으로 지반을 지탱했으나, 이후 점차 함몰 범위가 확대되었다(Gayle, 2025a; 2025b; Hashimoto, 2025; NHK, 2025). 따라서 Run4의 결과는 단기적인 지표 안정성만으로는 안전성을 평가할 수 없으며, 공동 내부의 잠재적 잔류 공간까지 고려한 장기적 대책이 필수적임을 시사한다.

Future work

본 연구는 MPM의 기본 구조를 활용하여 지하 공동 확대에 따른 붕괴 과정을 단순화된 지형 조건에서 모의한 것으로, 향후 과제에서는 다음과 같은 한계를 극복 개선하는 연구가 필요할 것으로 생각된다.

본 연구에서 제시된 수치모의는 단순화된 지반 구조와 경계 조건 하에서 수행되었으며, 실제 도심지의 복잡한 지반 조건과 매설물 환경을 충분히 반영하지 못하였다. 따라서 제시된 결과는 붕괴 메커니즘의 개념적 이해를 제공하는 수준으로 해석되어야 하며, 실규모 현장 상황에 대한 직접적인 예측에는 한계가 있다.

또한 본 모형은 현장 계측자료나 실내 모형실험 결과와의 정량적 검증 과정을 거치지 못하였다. 실제 사례와 정성적으로 일치하는 거동은 확인되었으나, 수치적 예측 정확성을 확보하기 위해서는 추가적인 검증 과정이 필수적이다.

이를 통해 본 연구에서 제시한 개념적 모의가 실제 도시 싱크홀 위험평가 및 예방적 관리 대책 수립에 보다 직접적으로 기여할 수 있을 것이다.

Conclusion

본 연구에서는 MPM을 활용하여 지하 공동의 확대에 따른 도심지 싱크홀 발생 메커니즘을 모의하였다. 수치모의 결과, 공동 직경이 증가할수록 지표면 침하 폭은 감소하였으나, 내부에서는 잔류 공동과 수직 변형이 지속되는 경향을 확인하였다. 특히 Run4 조건에서는 지표면 침하가 제한적으로 나타났음에도 불구하고 공동이 완전히 충전되지 않아 향후 공동 확대나 2차 침하 가능성이 상존함을 보여주었다. 이는 실제 후쿠오카(2016) 사고의 급격한 전면 붕괴와 사이타마(2025) 사고의 국부적 초기 함몰 사례와도 부합하는 결과이다. 따라서 단기적 보수만으로는 안전성을 확보하기 어렵고, 내부 공동의 잠재적 잔류와 장기적 거동을 함께 고려한 정밀 조사 및 보강 대책이 필요하다. 본 연구는 도시 기반시설의 지반 안전성 평가 및 싱크홀 예방 관리에 있어 MPM 기반 수치해석이 유효한 도구임을 입증하며, 향후 관로 노후화 지역의 사전 위험도 평가 및 정책적 의사결정에 기초 자료로 활용될 수 있을 것이다.