MODFLOW

1. 개요[편집]

MODFLOW Modular Ground-Water Flow Model
개발	미국 지질조사국(USGS)
이산화	블록중심 유한차분 (= 사각격자 유한체적)
미지수	셀 중심의 수두 h 하나
구조	패키지(package) 모듈 조립식
계보	88 → 96 → 2000 → 2005 → 6 (+ NWT, USG)
언어·라이선스	Fortran · 퍼블릭 도메인
생태계	ModelMuse, FloPy, PEST, MODPATH, MT3D

지하수 하는 사람들이 “모델 돌렸다”고 하면 되묻지 않는다. 어차피 MODFLOW다.

MODFLOW는 미국 지질조사국(USGS)이 개발한 3차원 지하수 유동 유한차분 모델이자, 그 입출력 형식이 사실상 업계 표준이 되어 버린 코드다. 1980년대에 나온 이래 지하수 수치모델링의 기본값 자리를 한 번도 내준 적이 없고, 규제·소송·인허가 문서에서 “MODFLOW로 해석했다”는 문장은 별도의 정당화를 요구받지 않는다.

전산유체역학 쪽 사람에게 설명하자면 지하수판 OpenFOAM 이 가장 가까운 비유다. 무료·공개·모듈형이고, 상용 GUI들이 전부 이 엔진을 감싸는 껍데기이며, 논문의 절반이 이 위에서 돌아간다. 다만 결정적인 차이가 둘 있다. 첫째, OpenFOAM이 이산화 도식부터 갈아 끼우는 툴박스라면 MODFLOW는 이산화가 하나로 고정된 모델이다. 사용자가 조립하는 것은 도식이 아니라 물리 현상(우물·하천·배수)이다. 둘째, USGS 소프트웨어는 GPL이 아니라 퍼블릭 도메인이라 상용 벤더가 그대로 가져다 팔 수 있다. 이 라이선스 한 줄이 ModelMuse부터 Groundwater Vistas까지의 생태계를 만들었다.

풀려는 방정식 자체는 지하수 유동 문서에 있는 그 확산식이다. MODFLOW의 정체성은 방정식이 아니라 그 방정식을 어떤 자료구조로 자르고, 무엇을 패키지로 뽑아냈으며, 자유면 셀이 마를 때 무슨 일이 벌어지는가에 있다.

2. 블록중심 유한차분[편집]

격자는 행·열·층으로 인덱스된 직육면체 셀이고, 미지수는 각 셀 중심의 수두 하나다. 셀 $m$ 에 대한 이산식은 인접 6개 셀과의 흐름 + 외부 응력 = 저류 변화라는 물수지 그 자체다.

$$\sum_{n\in N(m)}C_{mn}\left(h_n-h_m\right)+P_mh_m+Q_m =S_{s,m}\,\Delta V_m\,\frac{h_m^{t}-h_m^{t-\Delta t}}{\Delta t}$$

여기서 $C_{mn}$ 이 전도도(conductance)이고, 이것이 MODFLOW의 심장이다. 행 방향으로 이웃한 두 셀 사이는

$$CR_{i,j+1/2,k}=\frac{2\,\Delta c_i\,\Delta v_k}{\dfrac{\Delta r_j}{K_{i,j,k}}+\dfrac{\Delta r_{j+1}}{K_{i,j+1,k}}}$$

즉 셀 중심에서 면까지의 저항 두 개를 직렬로 더한 형태 — 수리전도도의 거리가중 조화평균이다. 산술평균이 아닌 이유는 수리전도도 문서의 $K_\perp$ 논의와 같다. 셀 경계를 가로지르는 흐름은 직렬 저항이고, 저투수 셀 하나가 흐름을 막는 물리를 산술평균은 재현하지 못한다.¹

여기서 중요한 사실 하나. 이 “유한차분”은 사실상 사각격자 위의 유한체적법이다. 셀 단위 물수지에서 출발했으므로 질량보존이 이산 수준에서 정확히 성립하고, 그래서 MODFLOW는 모의가 끝날 때마다 전 영역 물수지표를 찍어 불일치율(percent discrepancy)을 보고한다. 국룰은 1% 이하이며, 잔차는 얌전한데 이 값이 튀면 해가 아니라 모델 설정이 잘못된 것이다.

행렬은 대칭 양정치 희소행렬이 되고 — 나중에 볼 뉴턴 정식화만 아니면 — 대각지배가 보장된다. 수두의존 경계가 대각항에 양수를 더해 주기 때문에, 경계조건이 많을수록 오히려 잘 풀린다.

시간 축은 응력기간(stress period) → 시간스텝(time step) → 외부 반복 → 내부 솔버 반복의 4중 구조다. 응력기간은 양수량·함양량 같은 경계 조건이 일정하게 유지되는 구간이고, 그 안을 배수 $TSMULT$ 로 점점 키우며 쪼갠 것이 시간스텝이다. 시간적분은 완전 음해(후방오일러)만 쓴다 — 정확도를 포기하고 안정성을 산 선택이며, 40년 동안 아무도 이 결정을 되돌리지 않았다.

3. 패키지 구조[편집]

MODFLOW의 진짜 발명품은 수치도식이 아니라 패키지라는 모듈 개념이다. 물리 현상 하나마다 파일 하나·서브루틴 한 벌이 대응하고, 각 패키지가 하는 일은 위 이산식의 $P_m$(대각 기여)과 $Q_m$(우변 기여)을 채워 넣는 것뿐이다.

패키지	하는 일	이산식 기여
WEL	지정 유량 우물	$Q$ 만 (수두와 무관)
RCH	상부 함양	$Q$ 만
GHB	일반수두경계	$P$ 와 $Q$ 둘 다, 항상 활성
DRN	배수·용출	$h>z_{drain}$ 일 때만 활성
RIV	하천-대수층 교환	하상 아래로 떨어지면 유량 고정
EVT	증발산	수두 구간별 선형
CHD	시간에 따라 변하는 정수두	강제 디리클레
HFB	얇은 차수벽	면 전도도를 직접 깎음

여기서 GHB·DRN·RIV 셋은 전부 같은 제3종 경계조건 $q=C(H_{ext}-h)$ 이고, 차이는 스위치가 어떻게 붙느냐뿐이다. GHB는 양방향 상시, DRN은 배수고 위에서만 한쪽으로, RIV는 수두가 하상 바닥 아래로 내려가면 $q=C(H_{riv}-R_{bot})$ 로 고정 침투가 된다(하천과 대수층이 끊어져 물이 자유낙하하는 상황). 개념적으로는 같은 식이지만 DRN과 RIV는 해에 의존하는 스위치를 갖고 있어 비선형이며, 게다가 스위치 지점에서 미분이 꺾인다. 수렴이 안 될 때 제일 먼저 의심할 곳이 여기다.

내부 유동 물성 쪽은 계보를 따라 이름이 바뀌었다. BCF(블록중심 유동) → LPF(층 물성 유동) / HUF(수리지질단위) → UPW(상류가중, NWT 전용) → NPF(절점 물성, MODFLOW 6). 하는 일은 같다 — 층별 $K$, 이방성비, 피압/자유면 여부를 받아 전도도를 계산한다.

4. 마른 셀 — 이 문서의 킬포인트[편집]

MODFLOW를 실제로 돌려 본 사람이 가장 오래 붙잡고 있는 문제는 격자도 물성도 아니고 셀이 마르는 것(dry cell)이다.

자유면 셀의 투수량계수는 포화두께에 비례한다.

$$T=K\left(h-z_{bot}\right)$$

수두가 셀 바닥에 가까워지면 $T\to0$ 이고, 바닥 아래로 내려가면 물리적으로 의미가 없어진다. 고전적인 BCF/LPF 정식화는 여기서 셀을 죽인다 — IBOUND를 0으로 바꿔 행렬에서 행과 열을 통째로 빼고, 출력에는 사용자가 지정한 HDRY 값(관행적으로 $-10^{30}$)을 찍는다. 이 한 줄의 처리가 세 가지 재앙을 부른다.

첫째, 반복마다 푸는 문제가 달라진다. 픽카드 반복은 “같은 비선형 방정식의 고정점”을 찾는 알고리즘인데, 셀이 죽고 살아나면 미지수 집합 자체가 바뀐다. 수렴할 고정점이 존재하지 않는 상태에서 반복을 돌리는 것이므로, 잔차가 어떤 값 아래로 안 내려가고 두 상태 사이를 왕복한다. 전형적 증상: 최대 수두변화량이 몇백 반복째 똑같은 두 값을 번갈아 찍는다.

둘째, 마른 셀은 연직 흐름을 통째로 끊는다. 셀이 죽으면 그 위아래 연결도 사라진다. 상부 층의 함양이 아래 대수층으로 내려갈 통로가 막히므로, 양수 우물이 있는 하부 층이 공급원과 단절되어 수두가 더 떨어지고, 그 결과 이웃 셀이 또 마른다. 연쇄 건조가 모델 한 구석을 통째로 지워 버리는 일이 실제로 벌어진다. 최악은 이게 조용히 일어난다는 것 — 물수지는 맞고(없는 셀은 물도 안 쓰니까) 그림만 이상하다.

셋째, 재습윤(rewetting)은 문제를 더 키운다. WETDRY 옵션을 켜면 이웃 수두가 임계값을 넘을 때 마른 셀이 되살아나는데, 되살아난 셀에 넣을 수두는 추정치다. 그 추정이 너무 높으면 다음 반복에서 물이 빠져나가 다시 마르고, 낮으면 이웃이 밀어 넣어 다시 젖는다. 주기 2의 진동이 그대로 자리 잡고 수렴은 영원히 오지 않는다. 실무자들이 WETFCT, IWETIT 같은 계수를 손으로 만지작거리며 밤을 새우는 이유이며, 이 바닥에서 “수렴은 신에게 맡긴다”는 말이 나온 자리이기도 하다.²

근본 원인은 이산식이 $h$ 의 매끄러운 함수가 아니라는 것이다. 전도도가 어느 순간 계단식으로 0이 되고, 그 지점에서 미분이 존재하지 않는다. 뉴턴법을 쓰든 픽카드를 쓰든 매끄럽지 않은 함수의 영점을 찾는 일은 원리적으로 잘 안 된다.

처방은 2011년 MODFLOW-NWT와 UPW 패키지에서 정리됐고, MODFLOW 6의 NEWTON 옵션이 같은 아이디어를 이어받았다.

• 셀을 절대 죽이지 않는다. 대신 포화도 함수 $\Omega(h)$ 를 곱해 전도도를 0이 아닌 아주 작은 값까지 연속적으로 낮춘다. 미지수 집합이 고정되므로 문제가 반복 내내 동일해진다.
• $\Omega$ 를 매끄러운(미분 가능한) 3차 함수 같은 것으로 잡아 야코비안을 정의 가능하게 만든다.
• 면의 포화두께를 조화평균이 아니라 상류가중(upstream weighting)으로 잡는다. 다상유동에서 상대투과도를 상류가중하는 것과 같은 처방이고, 이것이 없으면 뉴턴 반복이 여전히 튄다.
• 뉴턴-랩슨법으로 외부 반복을 돌린다.

대가는 야코비안이 비대칭이 된다는 것이다. 대칭 양정치가 깨지므로 켤레기울기법을 못 쓰고 GMRES나 BiCGSTAB + ILU 계열로 갈아타야 하며, 반복당 비용이 오른다. 그래도 안 되던 문제가 되는 것과 빨리 안 되는 것 사이의 선택이라 아무도 망설이지 않는다. 지금은 자유면 층이 있으면 일단 뉴턴을 켜고 시작하는 것이 기본값에 가깝다.

5. 솔버와 격자의 계보[편집]

솔버. 초기 MODFLOW-88은 SIP(강한 음해 절차)와 SSOR를 썼다. 1990년대에 PCG2(켤레기울기법 + 다항식·불완전 촐레스키 전처리기)가 들어오면서 사실상 표준이 됐고, 이후 소규모 직접해법 DE4, 기하 다중격자법 기반 GMG, 비선형 처리를 강화한 PCGN이 추가됐다. MODFLOW-NWT는 앞서 본 이유로 비대칭 솔버(XMD, GMRES)를 물고 왔고, MODFLOW-USG의 SMS와 MODFLOW 6의 IMS가 이 흐름을 통합했다. 지금의 IMS는 비선형 반복(뉴턴 + 언더릴랙세이션 + 백트래킹)과 선형 반복을 한 블록에서 설정하는 구조다.

격자. 이것이 MODFLOW 6의 가장 큰 변화다.

• DIS — 전통적인 행·열·층 구조격자. 좌표만으로 이웃이 결정된다.
• DISV — 평면은 임의 다각형(보로노이 등), 연직은 층으로 쌓는 반구조 격자. 우물이나 하천 주변만 국소적으로 조밀하게 깔 수 있다.
• DISU — 완전 비구조. 셀과 연결(connection) 목록을 직접 주며, 수직 세분화(nested refinement)와 층 소멸(pinch-out)을 자연스럽게 표현한다. MODFLOW-USG(2013)가 도입한 제어체적 유한차분(CVFD) 정식화가 뿌리다.

비구조로 가면 대가가 따라온다. 2점 플럭스 근사는 셀 중심을 잇는 선분이 면과 직교할 때만 정합적이라, 뒤틀린 격자나 완전 텐서 이방성에서는 격자를 조여도 틀린 값으로 수렴한다. MODFLOW 6이 XT3D 옵션으로 다점 근사를 제공하는 것이 이 문제에 대한 답이다. 이 구조적 사정은 다르시 법칙 문서의 TPFA·MPFA 논의와 정확히 같다.

버전 계보를 한 줄로. 88(모듈 개념 확립) → 96(파일 구조 정리) → 2000(프로세스 개념 도입, 민감도·매개변수 추정 내장) → 2005(포트란 모듈 기반 메모리 관리로 재작성) → 6(2017, 객체지향 재작성 · 비구조 격자 · 하나의 모의 안에 여러 모델을 두고 교환(exchange)으로 잇는 구조). 곁가지로 NWT(2011, 자유면 뉴턴), USG(2013, 비구조), 밀도 연성 SEAWAT이 있고, MODFLOW 6에는 뒤이어 용질 수송 모델(GWT)과 열 수송 모델이 합류해 예전에 MT3DMS로 외주 주던 일을 같은 실행 파일 안에서 처리한다.

6. 생태계 — 코드보다 큰 것[편집]

MODFLOW의 지배력은 코드 품질보다 주변부에서 나온다. 텍스트 입력 파일 형식이 표준이 되면서, 그 형식을 읽고 쓰는 도구들이 하나의 산업을 이뤘다.

• ModelMuse — USGS가 만든 무료 GUI. 격자와 무관하게 도형(점·선·면)으로 물성과 경계를 정의하고 나중에 격자에 투영하는 방식이라, 격자를 갈아엎어도 입력을 다시 안 만들어도 된다.
• FloPy — 파이썬으로 입력을 생성하고 결과를 읽는 패키지. 지금 나오는 MODFLOW 논문의 재현 가능한 스크립트는 사실상 전부 이것이다. 수천 개 시나리오를 돌리는 앙상블 작업에서 GUI는 애초에 선택지가 아니다.
• PEST / PEST++ — 모델 독립적 보정·불확실성 도구. 관측 수두와 유량에 맞춰 물성을 조정하는 역문제를 풀며, 파일럿 포인트를 크리깅으로 보간하고 티코노프 정칙화를 거는 실무 표준 절차가 여기 들어 있다. 앙상블 계열(PESTPP-IES)은 앙상블 칼만 필터의 반복형 변종이다.
• MODPATH — 셀 면 유량에서 속도장을 복원해 입자를 추적한다. 우물의 포획구역(capture zone)을 그리는 표준 도구.
• 상용 GUI — Groundwater Vistas, GMS, Visual MODFLOW Flex. 퍼블릭 도메인 엔진 위에 지리정보·시각화·보정 워크플로를 얹어 판다.

7. 한계[편집]

• 격자가 물리를 제약한다. DIS 시절의 직교 격자는 하천을 계단으로 근사하고, 층이 얇아지는 지점을 표현하지 못했다. 6의 DISV/DISU가 이 문제를 상당히 풀었지만, 기존 자산의 대부분은 여전히 구조격자 위에 있다.
• 불포화대가 없다. MODFLOW는 포화 유동 코드다. UZF 패키지가 통기대를 간이 저류 모형으로 흉내 낼 뿐, 리처즈 방정식을 푸는 것이 아니다. 침투 지연이나 모세관대가 중요한 문제에는 애초에 적합한 도구가 아니다.
• 입력의 불확실성이 지배한다. 이건 코드 탓이 아니다. 수리전도도 문서가 말하는 그대로, 예측 오차는 이산화 오차보다 두세 자릿수 크다. 격자를 두 배로 조이는 노력이 시추공 하나만 못한 상황이 이 분야의 기본값이며, 그래서 검증 및 확인과 민감도 해석이 코드 성능보다 훨씬 큰 비중을 차지한다.
• 수렴 실패의 절반은 사용자 책임이다. 마른 셀, 켜졌다 꺼지는 배수 셀, 자릿수가 어긋난 전도도, 응력기간 경계의 계단 변화 — 원인 목록이 거의 정해져 있는데도 매번 새로 겪는다.³

8. 관련 문서[편집]

• 지하수 유동 · 다르시 법칙 · 수리전도도
• 유한차분법 · 유한체적법 · 혼합유한요소법 · 격자
• 뉴턴-랩슨법 · 켤레기울기법 · GMRES · 전처리기 · 다중격자법
• 희소행렬 · 경계 조건 · 대류-확산 방정식
• FloPy · MODPATH · MT3DMS
• 역문제 · 크리깅 · 앙상블 칼만 필터 · 민감도 해석
• OpenFOAM · EPANET · HEC-RAS · 검증 및 확인

9. Footnotes[편집]

1. 이 조화평균 하나가 초보자를 반드시 한 번은 넘어뜨린다. 저투수 층을 표현하려고 셀 하나에 $K$ 를 $10^{-9}$ 로 박아 놓으면, 그 셀과 이웃 사이의 전도도가 이웃 $K$ 와 무관하게 거의 0이 된다. 의도한 대로 물이 막히니 다행 같지만, 동시에 행렬의 조건수가 수직 상승해 솔버가 기어가기 시작한다. 차수층은 셀 물성으로 표현하는 대신 HFB 패키지로 면 전도도만 깎는 것이 정석인 이유. ↩

2. 재습윤 옵션을 켜는 것이 언제나 나쁘다는 뜻은 아니다. 다만 매뉴얼조차 사실상 “가급적 뉴턴 정식화를 쓰라”는 쪽으로 정리됐다. 2011년 이전에 작성된 모델을 물려받아 WETDRY 배열을 손보고 있다면, 그 시간에 UPW/NPF로 옮기는 것이 대체로 더 빠르다. 물론 인허가 문서에 제출한 모델의 버전을 바꾸려면 심의를 다시 받아야 한다는 함정이 있다. 그래서 2005 시절 모델이 아직도 돌아간다. ↩

3. 수렴이 안 될 때 폐쇄기준(HCLOSE, RCLOSE)을 슬쩍 느슨하게 푸는 유혹이 늘 있다. 그러면 실행은 끝나지만 물수지 불일치율이 같이 올라가고, 보고서 어딘가에 “불일치율 3.7%“가 조용히 박힌다. 심의위원이 그 표를 볼 확률은 낮지만, 그 표를 보는 심의위원은 반드시 그것부터 본다. ↩