지하수 유동

유체역학 전산유체역학 시뮬레이션 마지막 수정: 2026-08-29 04:14:52

1. 개요[편집]

지하수 유동 Groundwater flow
지배 방정식	S_s ∂h/∂t = ∇·(K∇h) + W
미지수	수두 h [m] — 압력이 아니다
2차원 형태	S ∂h/∂t = ∇·(T∇h) + W, T = Kb, S = S_sb
정상상태	K 균질이면 ∇²h = 0 (라플라스)
자유면	부시네스크 식 — T = Kh 라 비선형
대표 해석해	티엠(1906, 정상) · 타이스(1935, 비정상)
대표 코드	MODFLOW, FEFLOW

지표수는 어디로 흐르는지 눈으로 보인다. 지하수는 안 보이는데, 그래서 방정식이 필요하다.

지하수 유동은 포화된 다공성 매질 안에서 물이 수두 구배를 따라 이동하는 현상으로, 다르시 법칙에 질량보존을 결합해 얻은 확산형 편미분방정식으로 기술된다. 미지수는 압력이 아니라 수두(hydraulic head) $h = z + p/\rho g$ 다. 지하수는 압력이 높은 곳이 아니라 수두가 높은 곳에서 낮은 곳으로 흐르며, 이 구별을 놓치면 깊은 관정이 얕은 관정보다 압력이 높은데도 물이 위로 솟는 자분(自噴) 현상이 설명되지 않는다.¹

매질 물성 $\mathbf{K}$ 자체는 수리전도도 문서에, 다르시 법칙의 유도와 적용 한계는 다르시 법칙 문서에 있다. 이 문서는 지배식·경계조건·해석해·수치해법을 다룬다.

2. 지배 방정식[편집]

검사체적 하나에 대해 유입-유출 = 저류 변화를 쓰고 플럭스를 다르시 법칙으로 대체하면

S_s\frac{\partial h}{\partial t}=\nabla\cdot\left(\mathbf{K}\nabla h\right)+W

를 얻는다. $W$ 는 단위 부피당 원천·소멸항(양수, 함양, 누수)이고, 비저류계수 $S_s\,[\mathrm{m^{-1}}]$ 는 수두가 1 m 내려갈 때 단위 부피에서 빠져나오는 물의 부피다. 야콥의 분해로

S_s=\rho g\left(\alpha+n\beta\right)

$\alpha$ 는 골격의 압축률, $\beta$ 는 물의 압축률( $\approx4.4\times10^{-10}\,\mathrm{Pa^{-1}}$ ), $n$ 은 공극률이다. 물이 나오는 이유가 물이 팽창해서가 아니라 골격이 압축돼서라는 것이 이 식의 요점이며, 그래서 과잉 양수는 지반 침하로 이어진다. 형태만 보면 이 식은 비등방 열전도 방정식과 글자 그대로 같은 골격이고, 그 덕분에 20세기 초 지하수학자들이 열전도 문헌을 통째로 번역해 쓸 수 있었다.

실무에서는 대수층 두께 $b$ 로 연직 적분해 2차원으로 낮춘다. 이때 투수량계수 $T=Kb\,[\mathrm{m^2/s}]$ 와 저류계수 $S=S_sb$ (무차원)가 등장한다.

S\frac{\partial h}{\partial t}=\nabla\cdot\left(T\nabla h\right)+W

3. 피압과 자유면 — 그리고 비선형성[편집]

구분	피압대수층	자유면대수층
상부 경계	난투수층에 눌려 있음	지하수면(대기압)
저류 기구	골격·물의 압축	공극의 실제 배수
저류계수	$S \sim 10^{-5}\!-\!10^{-3}$	$S_y \sim 0.01\!-\!0.3$
투수량계수	$T=Kb$ , 상수	$T=Kh$ , 미지수의 함수
방정식	선형	비선형

피압대수층에서 $T$ 는 상수다. 수두가 내려가도 대수층 두께는 그대로이므로(물은 여전히 꽉 차 있고 압력만 낮아진다) 방정식이 깨끗한 선형 확산식이 된다.

자유면대수층은 사정이 다르다. 수두 자체가 포화대의 상단이므로, 바닥에서 잰 수두 $h$ 가 곧 포화두께이고 $T=Kh$ 가 된다. 두푸이 가정(흐름이 수평이고 연직으로 정수압 분포)을 얹으면 부시네스크 방정식이 나온다.

S_y\frac{\partial h}{\partial t}=\nabla\cdot\left(K\,h\,\nabla h\right)+N

$\nabla\cdot(Kh\nabla h)=\tfrac{K}{2}\nabla^2(h^2)$ 이므로 $h^2$ 를 미지수로 삼으면 정상상태에서는 선형이 된다. 1차원 정상 흐름에서 $h^2$ 가 위치에 대해 선형이라는 두푸이 포물선이 여기서 나온다.

q=\frac{K\left(h_1^2-h_2^2\right)}{2L}

여기에 재미있는 사실이 하나 붙는다. 두푸이 가정은 연직 유속을 통째로 버리므로 자유면의 모양은 틀리게 준다(하류 끝의 침윤면을 재현하지 못한다). 그런데 위 유량 공식만큼은 두푸이 가정 없이도 정확하다는 것이 알려져 있다. 수두 프로파일은 틀리는데 유량은 맞는 근사 — 실무가 이 가정을 60년 넘게 놓지 않는 이유다.

비정상 문제에서는 비선형이 그대로 남고, 여기에 저류계수 차이가 겹친다. $S_y$ 는 $S$ 보다 두세 자릿수 크므로 수리확산계수 $T/S$ 는 피압 쪽이 훨씬 크다. 피압대수층에서 $T=10^{-3}\,\mathrm{m^2/s}$ , $S=10^{-4}$ 이면 $T/S=10\,\mathrm{m^2/s}$ , 자유면에서 $S_y=0.2$ 면 $5\times10^{-3}\,\mathrm{m^2/s}$ — 압력 신호가 피압대수층에서 2000배 빨리 퍼진다. 같은 양을 뽑아도 피압대수층은 수 km 밖까지 얕게 수위를 떨어뜨리고, 자유면대수층은 관정 근처만 깊게 파낸다. 이웃 마을에서 민원이 오는 쪽은 언제나 전자다.

4. 정상상태와 유선망[편집]

$\partial h/\partial t=0$ 이고 원천항이 없으면 남는 것은

\nabla\cdot\left(\mathbf{K}\nabla h\right)=0 \quad\xrightarrow{\ \mathbf{K}=\text{const}\ }\quad \nabla^2h=0

즉 라플라스 방정식이다. 정상 지하수 유동은 정상 열전도, 정전기 퍼텐셜, 포텐셜 유동과 수학적으로 완전히 같은 문제이며, 조화함수의 최대값 원리·평균값 성질·유일성이 그대로 성립한다. 실무적으로는 정상해에는 $K$ 의 절대 크기가 들어가지 않는다는 사실이 중요하다. $\mathbf{K}\to c\mathbf{K}$ 로 통째로 키워도 같은 $h$ 가 나오므로, 수두만 맞춰서 보정한 모델은 $K$ 를 상수배까지밖에 결정하지 못한다.

2차원 등방·균질 매질에서는 수두 $h$ 와 유동함수 $\psi$ 가 켤레조화함수이므로 등수두선과 유선이 직교하고, 이 직교 곡선망을 손으로 그린 것이 유선망(flow net)이다. 곡선사각형이 “정사각형에 가깝게” 되도록 그리면 유량이 칸 세기만으로 나온다.

Q=K\,H\,\frac{n_f}{n_d}

$H$ 는 전수두차, $n_f$ 는 유로(flow channel) 개수, $n_d$ 는 등수두 구간 수다. 댐 하부 침투량이나 차수벽 효과를 모눈종이로 계산하던 시절의 기술인데, 지금도 수치해 결과의 자릿수를 검산하는 데 쓴다. 이방성 매질에서는 좌표를 $x'=x\sqrt{K_z/K_x}$ 로 눌러 등방으로 만든 뒤 유선망을 그리고 다시 되돌린다. 직교성은 변환된 공간에서만 성립한다는 것을 잊고 원래 좌표에서 직각으로 그리는 것이 이 분야 답안지의 단골 감점 요인.

5. 우물 문제[편집]

티엠 해(Thiem, 1906) — 정상상태 방사 흐름. 반경 $r$ 인 원통을 통과하는 유량이 어디서나 $Q$ 로 같아야 하므로 $Q=2\pi r T\,\mathrm{d}h/\mathrm{d}r$ , 적분하면

h_2-h_1=\frac{Q}{2\pi T}\ln\frac{r_2}{r_1}

주목할 점은 이것이 두 관측점 사이의 수두 차로만 쓰였다는 것이다. 무한 대수층에서 로그는 발산하므로 절대 수두를 주는 정상해는 존재하지 않는다. 실제로 정상상태가 성립하려면 유한한 거리에 하천이나 호수 같은 함양 경계가 있어야 하며, 교과서의 “영향반경 $R$ “은 물리량이 아니라 그 경계를 손으로 끼워 넣은 자리다.²

타이스 해(Theis, 1935) — 비정상. 무한 균질 피압대수층에 선원(line sink)을 놓고 푼 결과로, 무한 평판의 순간 선열원 해를 그대로 옮겨 온 것이다.

s(r,t)=\frac{Q}{4\pi T}W(u),\qquad W(u)=\int_u^\infty\frac{e^{-y}}{y}\,\mathrm{d}y,\qquad u=\frac{r^2S}{4Tt}

$s$ 는 수위강하, $W(u)$ 는 우물함수(= 지수적분 $E_1$ )다. $u$ 가 작아지는 후기(먼 시간, 가까운 거리)에는 $W(u)\approx-\gamma-\ln u$ 로 잘라 쓸 수 있고, 이것이 쿠퍼-야콥 근사(1946)다.

s\approx\frac{Q}{4\pi T}\ln\frac{2.25\,Tt}{r^2S}

$u<0.01$ 이면 오차가 1% 아래다. 이 형태가 실무를 지배하는 이유는 계산이 쉬워서가 아니라 $s$ 대 $\log t$ 가 직선이기 때문이다. 반로그 모눈에 찍어 직선을 긋고, 한 로그 사이클당 기울기 $\Delta s$ 에서 $T=2.303Q/(4\pi\Delta s)$ , $s=0$ 절편 $t_0$ 에서 $S=2.25Tt_0/r^2$ 를 읽는다. 양수시험 해석이 100년 가까이 자와 연필로 가능했던 근거이며, 지금도 자동 곡선맞춤 결과가 이 직선과 안 맞으면 무언가 잘못된 것이다.

누수 대수층(반투수층을 통해 위아래에서 물이 새어 들어오는 경우)은 한투시-야콥 해로 확장된다. 누수항이 붙으면 수위강하가 무한정 커지지 않고 유한한 값으로 수렴하며, 이때 나타나는 누수인자 $B=\sqrt{Tb'/K'}$ 가 반투수층 물성을 역산하는 통로가 된다.

6. 경계조건[편집]

지하수 문제의 난이도는 방정식이 아니라 경계조건에서 나온다.

제1종(디리클레), 정수두. $h=h_D$ . 해안선, 큰 호수, 수위가 고정된 하천. 무한한 물을 공급하는 경계라 모델이 필요로 하는 만큼 물을 뽑아 준다 — 경계를 관심 영역 가까이 두면 양수 영향이 통째로 흡수되어 사라진다.
제2종(노이만), 지정유량. $-\mathbf{K}\nabla h\cdot\mathbf{n}=q_N$ . 특수한 경우가 $q_N=0$ 인 불투수 경계이고, 지하수 분수령도 같은 조건으로 처리한다. 다만 분수령은 양수를 시작하면 움직이므로, 고정 무유량 경계로 잡아 두면 장기 예측이 조용히 틀린다.
제3종(로빈·코시), 수두의존 유량. $q_n=C\,(h_{ext}-h)$ . 하상 퇴적층을 통한 하천-대수층 교환, 일반수두경계, 반투수층 누수가 전부 이 형태다. 이산화하면 대각항에 양수 $C$ 를 더하는 꼴이라 행렬의 대각지배를 강화한다. 즉 수치적으로는 고마운 경계다.
단방향 경계. 배수구·용출은 $h>z_{drain}$ 일 때만 물을 빼고 그 아래에서는 아무 일도 하지 않는다. 조건이 해에 의존하므로 비선형이고, 게다가 미분불연속이다. 뉴턴 반복이 두 상태 사이를 왕복하며 수렴하지 않는 전형적인 원인.
자유면. 지하수면은 미지의 위치에 있는 이동 경계이며 조건이 두 개 걸린다. $p=0$ (즉 $h=z$ )이라는 기하 조건과, 함양·배수를 반영한 운동학적 조건이다. 미지수인 경계 위에 조건 두 개 — 자유면 문제가 비선형인 근본 이유다. 여기에 하류 사면에서는 $h=z$ 이면서 유출만 허용되는 침윤면이 붙는데, 그 범위 자체가 해의 일부인 상보성 조건이다.

7. 수치해법의 관점[편집]

이산화. 사각 격자 위의 블록중심 유한차분이 사실상 표준이고(MODFLOW), 그 이유는 셀 사이 전도도를 조화평균으로 잡으면 유한차분이 곧 유한체적법과 같아져 셀 단위 질량보존이 정확히 성립하기 때문이다. 복잡한 경계나 강한 이방성에는 유한요소법을 쓰지만, 표준 갤러킨은 요소 경계에서 플럭스가 불연속이라 오염물 이송 솔버에 넘길 속도장으로 곧바로 못 쓴다. 그래서 혼합유한요소법이나 후처리 플럭스 복원이 붙는다.
시간. 완전 음해(후방오일러)가 기본값이다. 양해법을 쓰면 $\Delta t\le S\Delta x^2/(4T)$ 인데, 피압대수층 값 $T/S=10\,\mathrm{m^2/s}$ 와 $\Delta x=100\,\mathrm{m}$ 를 넣으면 $\Delta t\lesssim250\,\mathrm{s}$ — 50년짜리 모의를 돌리려면 600만 스텝이 필요하다. 크랭크-니콜슨법은 정확도는 높지만 양수 개시처럼 계단식 응력 변화에서 진동하므로 실무 코드는 대체로 후방오일러를 고수한다.
선형계. 얻어지는 행렬은 대칭 양정치 희소행렬이라 켤레기울기법 + 불완전 촐레스키가 국룰이고, 층 방향으로 격자 간격이 극단적으로 눌린 문제(수평 1 km, 연직 1 m)에서는 다중격자법이 유리하다.
비선형. 자유면과 수두의존 경계 때문에 매 시간스텝에 외부 반복이 붙는다. 픽카드 반복은 안정적이지만 느리고, 뉴턴-랩슨법은 빠르지만 위에서 본 미분불연속 때문에 매끄러움 처리를 해 줘야 한다. 이 대목이 MODFLOW 계보의 절반을 설명한다.
검산. 모의가 끝나면 전 영역 물수지를 뽑아 유입-유출-저류변화의 불일치율이 1% 아래인지 본다. 잔차가 내려갔다고 답이 맞는 것이 아니다 — 마른 셀이 생기거나 경계가 잘못 켜지면 잔차는 얌전한데 물수지가 깨진다.

8. 무엇이 진짜 어려운가[편집]

방정식은 선형 확산식이고 솔버는 40년째 성숙해 있다. 남는 난이도는 전부 입력에 있다.

$\mathbf{K}$ 장의 추정은 대놓고 불량조건인 역문제이며(수리전도도 참조), 공간 분포는 크리깅이나 조건부 시뮬레이션으로 채운다. 함양량은 직접 못 재고 강우·증발산·토지이용에서 추정한다. 대수층 바닥의 형상은 시추공 몇 개를 이은 추정이다. 그 결과 지하수 모델의 예측 불확실성은 수치오차보다 두세 자릿수 크고, 이 분야에서 불확실성 정량화와 앙상블 예측이 유난히 일찍 자리 잡은 이유가 된다.³

응용 쪽도 대개 유동 자체가 목적이 아니다. 오염물 거동은 여기서 얻은 속도장을 대류-확산 방정식에 넘겨 푸는 것이고, 해수 침투는 밀도가 수두에 되먹임되는 연성 문제, 지반 침하는 유효응력을 매개로 한 수리-역학 연성 문제다. 하천·관망 쪽 형제 분야로는 HEC-RAS·SWMM·EPANET이 있다.⁴

9. 관련 문서[편집]

압력수두만 보고 흐름 방향을 판단하는 실수는 연직 방향에서 반드시 터진다. 깊이 100 m 관정의 압력은 지표 근처보다 당연히 높지만, 그 물이 위로 올라오느냐 마느냐는 $z+p/\rho g$ 를 비교해야 알 수 있다. 두 관정의 수위가 같아도 심도가 다르면 그 사이에 연직 흐름이 없다는 뜻이지, 물이 안 흐른다는 뜻이 아니다. ↩
“영향반경”을 마치 대수층의 물성인 양 표에서 찾아 넣는 관행이 아직도 있다. 지크하르트 공식 $R=3000s\sqrt{K}$ 같은 것들인데, 출처를 거슬러 올라가면 특정 현장의 회귀식이다. 무한 대수층에서 정상해가 없다는 사실을 알고 나면 이 숫자의 정체가 “적당히 먼 곳”이라는 것이 보인다. 값을 두 배로 바꿔도 로그 안에 들어가므로 결과가 별로 안 변한다는 것이 그나마의 위안이자, 동시에 이 숫자에 신경 쓸 필요가 없다는 증거다. ↩
그래서 이 바닥에서는 “모델을 믿느냐”보다 “무엇을 예측하는 데 쓰느냐”를 먼저 묻는다. 상대적 비교(관정을 여기 뚫는 게 저기보다 나은가)는 꽤 잘 맞고, 절대 예측(2050년 수위는 몇 m인가)은 대체로 못 맞힌다. 보고서의 결론 문장이 늘 “~에 비해 ~할 것으로 예상된다”인 데는 이유가 있다. ↩
지하수 모델러가 지표수 모델러와 회의를 하면 시간 스케일에서 먼저 부딪힌다. 한쪽은 홍수 첨두를 분 단위로 보고, 다른 쪽은 응력기간을 월 단위로 쪼갠다. 지표수-지하수 통합 모델이 어려운 이유의 절반은 물리가 아니라 이 시간 규모 차이이고, 나머지 절반은 하상 전도도라는 아무도 측정한 적 없는 계수다. ↩