해수 침투

1. 개요[편집]

해수 침투 Seawater intrusion
대상	해안 대수층의 담수-염수 경계 거동
1차 근사	가이벤-헤르츠베르크 관계 z ≈ 40 h
정식화	급격경계면 vs 변수밀도 유동-이송 완전결합
결합 고리	ρ(C) 가 유동을, 유동이 C 를 결정 — 양방향
표준 벤치마크	헨리 문제(준해석해) · 엘더 문제(격자 의존)
대표 코드	SEAWAT · SUTRA · FEFLOW · MODFLOW 6 BUY
실무 함정	염분 1%만 넘어도 식수·관개가 끝난다 — 오차 허용폭이 작다

담수를 1 m 더 뽑으면 염수 경계가 40 m 올라온다. 이 40배가 해안 대수층 관리의 전부다.

해수 침투는 해안 대수층에서 밀도가 큰 염수가 담수 영역으로 내륙 방향 또는 상향으로 밀고 들어오는 현상이며, 수치해석의 관점에서는 변수밀도 지하수 유동과 용질 이송이 양방향으로 완전 결합된 비선형 문제다. 지하수 유동 문서가 다루는 표준 확산방정식은 밀도가 상수라는 가정 위에 서 있는데, 여기서는 그 가정이 깨진다. 밀도가 농도를 통해 유동을 바꾸고, 그 유동이 다시 농도를 실어 나른다.

물리적 어려움과 별개로 이 문제가 사회적으로 중요한 이유는 허용 오차가 유난히 작기 때문이다. 해수의 염분은 약 35 g/L, 식수 기준은 대체로 0.25~0.5 g/L 수준이므로, 해수가 1% 남짓만 섞여도 관정 하나가 폐공이 된다. 오염 농도 예측을 두 배 틀려도 되는 문제가 아니라는 뜻이고, 그 결과 이 분야는 유난히 벤치마크 문화가 발달했다.

2. 가이벤-헤르츠베르크 관계[편집]

가장 단순한 모형은 담수와 염수를 섞이지 않는 두 유체로 보고 그 사이에 급격 경계면(sharp interface)이 있다고 가정하는 것이다. 경계면이 정지해 있고 염수 쪽이 정체(靜止)해 있다고 하면, 경계면 위 한 점에서 양쪽 정수압이 같아야 한다. 해수면 기준으로 담수 수위가 $h$ 만큼 높고 경계면이 $z$ 만큼 깊다면

$$\rho_f\,g\,(h+z)=\rho_s\,g\,z \quad\Longrightarrow\quad z=\frac{\rho_f}{\rho_s-\rho_f}\,h$$

이다. $\rho_f=1000\,\mathrm{kg/m^3}$, $\rho_s=1025\,\mathrm{kg/m^3}$ 을 넣으면 $1000/25=40$ 이므로

$$z\approx40\,h$$

가이벤-헤르츠베르크 관계(Ghyben–Herzberg)다.¹ 담수 수위가 해수면보다 1 m 높으면 그 아래 40 m까지가 담수라는 뜻이고, 뒤집으면 양수로 수위를 1 m 낮출 때마다 경계면이 40 m 올라온다는 뜻이다. 해안 관정이 몇 년을 멀쩡히 돌다가 어느 해에 갑자기 짜지는 것이 이 지렛대 비 때문이다. 밀도차가 2.5%밖에 안 되는데 기하학적 증폭은 40배 — 이 분야에서 밀도차를 무시할 수 없는 이유가 한 줄로 설명된다.

한계도 분명하다.

• 연직 흐름을 무시한다(정수압 가정). 담수가 실제로 바다로 흘러나가는 해안선 근처에서 가장 크게 틀린다.
• 유출면을 못 만든다. 해안선에서 $h=0$ 이면 $z=0$ 이 되어 경계면이 해안선에서 지표와 만나야 하는데, 그러면 담수가 나갈 통로가 없다. 실제로는 해안선에서 조금 떨어진 해저에 유한한 폭의 담수 유출면이 있다. 글로버(Glover, 1959)의 해석해가 이 유출면 폭을 담수 유출량과 함께 주어 이 모순을 해결한다.
• 경계면이 없다. 실제로는 분자확산과 수리분산 때문에 두께 수 m~수십 m의 전이대(transition zone)가 존재하고, 조석과 계절 변동이 이를 두껍게 만든다. 급격경계면 모형은 전이대의 중앙(농도 50% 등농도선)을 근사하는 것으로 이해해야 한다.

그럼에도 이 관계가 살아남은 것은 관측 없이도 자릿수를 준다는 점 때문이다. 급격경계면 가정 위에서 스트랙(Strack, 1976)이 담수-염수 두 영역을 하나의 퍼텐셜로 묶는 정식화를 내놓아, 경계면 선단(toe)의 위치를 해석적으로 계산할 수 있게 했다. 지금도 스크리닝 도구로 쓴다.

3. 변수밀도 유동과 이송의 결합[편집]

제대로 풀려면 두 방정식을 같이 푼다. 유동 쪽은 다르시 법칙에 밀도를 되살린 형태이고, 실무 코드는 담수 등가수두 $h_f=p/(\rho_f g)+z$ 를 미지수로 쓴다.

$$\mathbf{q}=-\frac{k\rho g}{\mu}\left(\nabla h_f+\frac{\rho-\rho_f}{\rho_f}\nabla z\right)$$

괄호 안 두 번째 항이 부력항이며, 밀도가 균일하면 사라진다. 이송 쪽은 대류-확산 방정식이다.

$$\frac{\partial(\phi C)}{\partial t}+\nabla\!\cdot\!(\mathbf{q}C)=\nabla\!\cdot\!\left(\phi\,\mathbf{D}\nabla C\right)+q_sC_s$$

닫는 것은 상태식 $\rho=\rho_f(1+\epsilon C)$ 다. 해수 농도 35 kg/m³ 에서 밀도가 1000에서 1025로 가므로 $\partial\rho/\partial C=25/35\approx0.71$ 이고, 이 값이 SEAWAT 계열의 기본 밀도-농도 기울기로 그대로 들어간다.

결합의 성격이 중요하다. 밀도가 유동에 되먹임되므로, 이송을 한 번 풀 때마다 유동을 다시 풀어야 한다. 표준 전략은 시간스텝마다 유동-이송을 번갈아 반복하는 픽카드 반복(SEAWAT의 결합 반복)이고, 강한 밀도 대비에서는 이 반복 자체가 진동해 완화계수가 필요하다. 밀도-농도 관계까지 야코비안에 넣는 완전 뉴턴-랩슨법은 FEFLOW·SUTRA 계열 일부에서 지원하는데, 수렴은 빠르지만 야코비안 조립 비용이 만만찮다.

수치적 난점은 이송 쪽이 이류 지배라는 것이다. 격자 페클레 수 $\mathrm{Pe}=\Delta x/\alpha_L$ 가 2를 넘으면 중심차분이 진동하고, 진동으로 농도가 음수가 되면 밀도가 담수보다 가벼워져 없던 부력이 생긴다. 유동-이송이 결합돼 있기 때문에 이송 쪽의 사소한 오버슛이 곧바로 가짜 대류로 증폭된다는 것이 단일 이송 문제와의 결정적 차이다. 그래서 TVD 제한자, 특성선 기반(MOC/MMOC) 기법, 또는 눈물을 머금은 격자 세분이 필요하다.

또 하나, 보스와 소자(Voss & Souza, 1987)가 지적한 정합 속도 근사(consistent velocity approximation) 문제가 있다. 압력 구배항과 부력항을 서로 다른 차수로 이산화하면, 물리적으로 정지해 있어야 할 성층 상태에서 두 항이 정확히 상쇄되지 않아 가짜 순환류가 생긴다. 유사 이송의 웰밸런스드 문제, 자유표면 유동의 정수압 균형 문제와 정확히 같은 병이다 — 원천항과 플럭스가 서로를 지워야 하는데 이산화가 그걸 모른다.

4. 헨리 문제와 엘더 문제[편집]

이 분야의 벤치마크 두 개는 성격이 정반대이고, 그래서 둘 다 필요하다.

헨리 문제(Henry, 1964). 2 m × 1 m 직사각 단면 대수층에 왼쪽에서 담수를 일정 플럭스로 넣고 오른쪽을 해수 수두로 잡아 정상상태 염수 쐐기를 구한다. 확산계수를 크게(고정 상수로) 잡아 놓은 덕분에 푸리에-갈러킨 급수로 준해석해가 존재한다는 것이 이 문제의 전부다. 변수밀도 코드를 새로 짜면 제일 먼저 돌려 보는 문제.

그런데 이 문제에는 유명한 반전이 있다. 심프슨과 클레멘트(Simpson & Clement, 2004)가 헨리 문제를 밀도 결합을 켠 경우와 끈 경우로 나눠 풀어 보니 두 결과가 거의 같았다. 즉 표준 헨리 문제는 담수 플럭스가 커서 밀도 효과에 별로 민감하지 않고, 따라서 “변수밀도 코드의 검증”이라는 원래 목적에 부적합하다. 그래서 담수 유입을 절반으로 줄인 수정 헨리 문제가 제안되었고, 지금은 이쪽이 더 정직한 검증으로 취급된다. 벤치마크를 통과했다는 것이 코드가 옳다는 뜻이 아니라는 교훈의 표본 사례다.²

엘더 문제(Elder, 1967의 열대류 실험 → Voss & Souza, 1987의 용질 버전). 위쪽 절반에 염수원을 놓고 아래로 손가락(finger) 모양 밀도류가 내려가는 불안정 문제다. 문제의 파라미터는 레일리 수 하나로 요약되며, 다공성 매질 자연대류의 임계값이 $\mathrm{Ra}_c=4\pi^2\approx39.5$ 인 데 반해 엘더 문제는 400 수준이라 대류가 강하게 발달한다.

여기서 이 문제의 악명이 시작된다. 격자를 세분할수록 손가락 개수가 바뀌고, 중앙에서 흐름이 올라가는 해와 내려가는 해가 격자에 따라 갈린다. 처음에는 코드 버그로 의심됐지만, 요한센(Johannsen, 2003)이 분기 해석을 돌려 이 파라미터 영역에 여러 개의 안정 정상해가 실제로 공존한다는 것을 보였다. 즉 격자를 세분하는 행위는 오차를 줄이는 것이 아니라 어느 해 가지로 끌려갈지를 바꾸는 것이고, 하나의 답으로 수렴할 것을 기대하는 것 자체가 틀렸다.

이것이 왜 중요한가. 우리가 격자 수렴 지수로 오차를 추정할 때 암묵적으로 깔고 있는 가정은 “격자를 줄이면 유일한 해로 간다”는 것이다. 엘더 문제는 그 가정이 성립하지 않는 물리계가 실재함을 보여주는 표준 반례다. 밀도 불안정이 낀 문제에서 격자 수렴 그래프를 자랑스럽게 첨부하기 전에, 그 곡선이 수렴이 아니라 분기 이론의 가지 선택일 가능성을 한 번은 의심해야 한다. 레일리-테일러 불안정과 새프먼-테일러 불안정에서 똑같은 이야기가 반복되는 것도 우연이 아니다.

이 두 문제 외에 HYDROCOIN 암염돔 문제(Case 5)와 SEAWAT 매뉴얼의 회전 밀도류 문제 등이 관례적으로 함께 돌려진다.

5. 코드[편집]

코드	기관·성격	특징
SEAWAT	USGS, 무료	MODFLOW + MT3DMS 를 밀도로 묶은 것. V4에서 열·점성까지
SUTRA	USGS, 무료	유한요소법 기반, 포화-불포화·용질/열 겸용
FEFLOW	DHI, 상용	3D FEM, GUI가 강하고 지열까지
MODFLOW 6	USGS, 무료	GWT 모형 + BUY(부력)·VSC(점성) 패키지로 SEAWAT 기능을 흡수
OpenGeoSys · Codesa-3D	학계 오픈소스	연구용, 연성 물리 확장

SEAWAT의 설계 철학은 볼 만하다. 새 코드를 짜는 대신 이미 검증된 유동 코드와 이송 코드를 밀도 결합 루프로 묶었다. 그 덕에 MODFLOW 생태계의 전처리·후처리·패키지를 통째로 물려받았지만, 동시에 MODFLOW의 유한차분 격자 제약(마른 셀, 층 구조)도 통째로 물려받았다. MODFLOW 6가 GWF-GWT 모형을 한 프로그램 안에 두고 BUY 패키지로 결합을 처리하는 쪽으로 간 것은 이 이중 상속을 정리하려는 시도로 읽힌다.

6. 실무 — 쐐기는 어떻게 전진하고 어떻게 막는가[편집]

과잉 양수가 만드는 일은 두 가지다.

• 염수 쐐기 전진(lateral intrusion). 내륙 수두가 낮아지면 경계면 선단이 육지 쪽으로 이동한다. 시간규모가 수년~수십 년이라 눈치챘을 때는 이미 늦다. 그리고 되돌리는 데는 들어온 시간의 몇 배가 걸린다 — 확산으로 두꺼워진 전이대는 수두를 회복해도 빨리 안 씻긴다.
• 염수 상승원추(upconing). 관정 바로 아래에서 경계면이 국소적으로 솟는다. 여기에는 유용한 경험칙이 있다. 경계면의 상승량이 관정 스크린과 초기 경계면 사이 거리의 약 3분의 1을 넘으면 불안정해져 급격히 관정까지 치솟는다(슈모락-메르카도). 즉 안전 양수량은 “짜지기 직전까지”가 아니라 그보다 훨씬 보수적으로 잡아야 한다.³

대책은 넷으로 정리된다. (가) 양수 재배치·감축 — 가장 싸고 가장 정치적으로 어렵다. (나) 주입 장벽(injection barrier) — 해안선을 따라 주입정을 늘어놓아 수두 능선을 세운다. 로스앤젤레스 카운티의 웨스트코스트 분지 장벽이 1950년대부터 돌아가는 대표 사례. (다) 추출 장벽 — 반대로 염수를 뽑아 내 쐐기를 끌어당긴다. 처분 문제가 남는다. (라) 물리적 차수벽 — 확실하지만 비싸다. 실무에서는 이들을 섞고, 그 조합의 최적화를 역문제·최적화 루프로 푸는 것이 최근 흐름이다.

기후변화와의 연결에서는 흔한 오해를 하나 정리해야 한다. “해수면이 오르면 침투가 그만큼 심해진다”는 직관은 반만 맞다. 워너와 시먼스(Werner & Simmons, 2009)의 구분에 따르면, 내륙 수두가 고정된 수두 지배계에서는 해수면 상승이 곧바로 큰 쐐기 전진으로 이어지지만, 내륙에서 들어오는 담수 플럭스가 고정된 계에서는 해수면이 올라가도 담수 유출량이 그대로라 쐐기 전진이 놀랄 만큼 제한적이다. 실제 대수층이 어느 쪽에 가까운가는 지형·함양 조건이 결정한다. 그리고 두 경우 모두에서, 해수면 상승보다 함양 감소와 양수 증가가 훨씬 큰 요인이라는 것이 반복적으로 보고된다. 사람이 뽑는 양이 바다보다 무섭다.

7. 무엇이 진짜 불확실한가[편집]

방정식은 정해져 있고 코드는 성숙했다. 남는 것은 지하수 유동에서와 같은 이야기다.

• 종·횡분산지수 $\alpha_L,\alpha_T$. 전이대 두께를 사실상 혼자 결정하는데, 척도 의존성이 심해 실험실 값과 현장 값이 두 자릿수 다르다. 모델의 전이대 두께가 관측과 안 맞을 때 제일 먼저 손대는 손잡이이자, 그래서 물리량이 아니라 보정 파라미터로 전락한 대표 사례다.
• 비균질성. 점토 렌즈 하나가 국소 침투 경로를 통째로 바꾼다. 크리깅이나 다지점 지구통계로 채워 넣고 앙상블을 돌리는 것 말고는 방법이 없다.
• 경계조건. 해저 유출면, 조석 변동, 해안 하천과의 교환이 전부 시간에 따라 변하고 관측이 어렵다.
• 관측 자료. 다심도 관측정(multi-level well)이 없으면 전이대의 연직 구조를 아예 못 본다. 단일 스크린 관정의 전기전도도 하나로 3차원 경계면을 보정하겠다는 것은 역문제로서 대놓고 불량조건이며, 그래서 이 분야는 불확실성 정량화를 일찍 받아들였다.

8. 관련 문서[편집]

• 지하수 유동 · 다르시 법칙 · 수리전도도 · 다공성 매질 유동
• MODFLOW · SEAWAT · SUTRA · MT3DMS
• 대류-확산 방정식 · 페클레 수 · 수치분산 · 뉴턴-랩슨법
• 자연대류 · 레일리-테일러 불안정 · 새프먼-테일러 불안정 · 분기 이론
• 유한요소법 · 유한차분법 · 혼합유한요소법 · 격자 수렴 지수
• 역문제 · 크리깅 · 불확실성 정량화 · 검증 및 확인
• 리처즈 방정식 · 유사 이송

9. Footnotes[편집]

1. 네덜란드의 바둔 가이벤(1888)과 독일의 알렉산더 헤르츠베르크(1901)가 독립적으로 발견했다. 두 사람 다 수리학자가 아니라 각각 육군 공병 장교와 상수도 기술자였고, 논문이 아니라 실무 보고서에 실었다. 정수압 두 줄이면 유도되는 관계를 두 나라가 13년 차이로 따로 발견했다는 것이 이 분야가 얼마나 현장에서 시작했는지를 보여준다. 참고로 밀도를 1000/1030으로 잡으면 계수가 33이 되므로, “40”은 물성이 아니라 관행값이다. ↩

2. 이 사건 이후 벤치마크 논문의 서술 방식이 조금 바뀌었다. “우리 코드가 헨리 해와 일치한다”만 쓰는 대신 “밀도 결합을 껐을 때와 켰을 때 결과가 얼마나 다른가”를 함께 보고하는 관행이 생겼다. 검증이란 결국 틀렸을 때 틀렸다고 말해 주는 시험이어야 한다는 당연한 이야기인데, 그걸 40년 걸려 깨달았다. ↩

3. 이 3분의 1 규칙이 실무에서 특히 고약한 이유는, 그 전까지는 아무 징후도 없다는 데 있다. 관정 수질이 몇 년간 완벽하게 깨끗하다가 어느 갈수기에 갑자기 전기전도도가 치솟고, 양수를 멈춰도 안 돌아온다. 상승원추는 비가역적이지 않지만 회복 시간규모가 형성 시간규모보다 훨씬 길다. “지금까지 괜찮았으니 앞으로도 괜찮다”가 통하지 않는 대표적인 지하수 문제이며, 그래서 해안 관정은 심도를 얕게 잡고 여러 개로 분산하는 것이 정석이다. ↩