SUTRA

소프트웨어 시뮬레이션 유체역학 마지막 수정: 2026-09-14 04:14:22

1. 개요[편집]

SUTRA Saturated-Unsaturated TRAnsport
개발	미국 지질조사국(USGS), 클리퍼드 보스(C. I. Voss) 외
초판	1984
이산화	혼합 유한요소 + 적분유한차분
요소	2차원 사각형 · 3차원 육면체
미지수	압력 p + 용질농도 또는 온도 U
포화 상태	포화 · 불포화 모두
결합	유동-수송 순차 반복(한 실행 파일 안에서)
언어·라이선스	Fortran · 퍼블릭 도메인

MODFLOW 진영이 “검증된 코드 둘을 묶자”고 할 때, SUTRA는 “처음부터 하나로 짜면 될 것을”이라고 답했다. 40년째 둘 다 살아 있다.

SUTRA는 USGS가 배포하는 밀도 의존 변포화 지하수 유동 및 용질·열 수송 해석 코드로, 유동과 수송을 서로 다른 두 코드에 외주 주지 않고 하나의 유한요소 코드 안에서 순차 반복으로 푸는 것이 특징이다. 이름 자체가 Saturated-Unsaturated Transport의 약자이고, 불포화대를 정직하게 푼다는 그 한 글자가 SEAWAT과의 결정적 차이다.

MODFLOW 계열이 사각격자 유한차분·패키지 조립식·수두 기반이라면, SUTRA는 유한요소·단일 코드·압력 기반이다. 설계 철학이 정반대이고, 그래서 두 코드를 비교하는 것은 밀도 결합 지하수 모델링의 선택지를 이해하는 가장 빠른 길이다. 풀려는 물리 자체는 해수 침투와 다공성 매질 유동 문서에 있고, 이 문서는 SUTRA가 그 물리를 어떤 정식화로 이산화하는가만 다룬다.

2. 두 개의 방정식[편집]

SUTRA가 푸는 것은 딱 두 개다. 유체 질량 보존식:

\left(S_w\rho S_{op}+\varepsilon\rho\frac{\partial S_w}{\partial p}\right)\frac{\partial p}{\partial t} +\varepsilon S_w\frac{\partial\rho}{\partial U}\frac{\partial U}{\partial t} =\nabla\!\cdot\!\left[\frac{\mathbf{k}\,k_r\,\rho}{\mu}\left(\nabla p-\rho\mathbf{g}\right)\right]+Q_p

그리고 통합 수송식:

\frac{\partial\left(\varepsilon S_w\rho U\right)}{\partial t} +\left(1-\varepsilon\right)\rho_s\frac{\partial\left(c_s U\right)}{\partial t} =\nabla\!\cdot\!\left[\varepsilon S_w\rho\left(\mathbf{D}_m+\mathbf{D}_d\right)\nabla U\right] -\nabla\!\cdot\!\left(\varepsilon S_w\rho\,\mathbf{v}\,U\right)+Q_U

$\varepsilon$ 는 공극률, $S_w$ 는 포화도, $k_r$ 는 상대투과도, $S_{op}$ 는 비저류율, $U$ 는 용질농도이거나 온도다.

여기서 두 가지를 짚어야 한다.

첫째, 미지수가 수두가 아니라 압력이다. 밀도가 변하면 수두는 애초에 잘 정의되지 않는 양이라, SEAWAT은 등가담수수두라는 환산 변수를 발명해서 MODFLOW의 자료구조를 재활용했다. SUTRA는 그 우회를 하지 않고 압력을 그대로 미지수로 삼는다. 대가는 관측 자료를 넣을 때 수위를 압력으로 환산해야 한다는 것이고, 이득은 밀도가 변해도 방정식이 변형되지 않는다는 것이다. 부력은 $\nabla p-\rho\mathbf{g}$ 안에 이미 들어 있다.

둘째, 수송식이 하나다. $U$ 를 농도로 읽으면 용질 수송이고 온도로 읽으면 열 수송이다. 열의 경우 고체 골격의 비열이 좌변 둘째 항(= 용질에서는 선형 흡착에 의한 지연)으로, 열전도가 $\mathbf{D}_m$ (= 용질에서는 분자확산)으로 들어간다. SEAWAT V4가 “온도를 화학종 하나로 취급”해서 열을 넣은 것과 수학적으로 같은 사상인데, SEAWAT은 그걸 우회로로 썼고 SUTRA는 처음부터 정식화에 넣어 두었다. 대신 SUTRA는 한 번에 하나만 푼다 — 용질과 열을 동시에 풀려면 다중종 파생판(SUTRA-MS)이 필요하다.¹

3. 혼합 유한요소-적분유한차분[편집]

SUTRA의 이산화는 순수한 갤러킨 방법이 아니다. 매뉴얼이 스스로를 “hybrid finite-element and integrated-finite-difference method”라고 부르는데, 이 어정쩡한 조합에는 이유가 있다.

플럭스 항(다르시 흐름, 분산, 이류)은 유한요소법으로 다룬다. 사각형·육면체 요소 위의 쌍선형/삼선형 기저함수와 갤러킨 가중잔차. 비스듬한 층리와 곡면 경계를 격자 계단으로 근사하지 않아도 되는 것이 유한요소를 택한 이유다.
저류 항과 원천 항은 절점별로 집중(lumping)시킨다. 즉 그 항만 적분유한차분처럼 다룬다. 일관 질량행렬을 그대로 쓰면 급격한 전선 근처에서 비물리적 진동이 생기고, 불포화 문제에서는 그 진동이 $S_w$ 를 1 위나 0 아래로 밀어내 계산이 죽는다. 질량 집중은 정확도를 조금 내주고 단조성을 사는 거래이고, 변포화 코드에서는 거의 예외 없이 채택된다.

이류항에는 상류가중 옵션이 있지만 매뉴얼이 대놓고 만류한다. 이유는 수치분산 문서 그대로 — 상류가중이 만드는 인위적 분산이 물리적 분산지수를 압도해 버리면 계산한 전이대 두께가 물성이 아니라 격자의 함수가 된다. SUTRA의 권고는 격자 페클레 수를 4 이하로 잡고 상류가중을 끄는 쪽이다.

4. 정합 속도 — 유한요소가 밀도류에서 한 번은 넘어지는 곳[편집]

SUTRA 문헌에서 가장 자주 인용되는 기술적 기여는 정합 속도(consistent velocity) 계산이다. 보스와 소자(Voss & Souza, 1987)가 정리한 이 문제는 이렇다.

정수압 평형 상태의 성층 밀도 분포를 생각하자. 물리적으로는 $\nabla p=\rho\mathbf{g}$ 가 정확히 성립해서 속도가 0이어야 한다. 그런데 요소 안에서 압력을 쌍선형 기저로 보간하면 $\nabla p$ 는 각 방향으로 차수가 하나 낮아진 함수가 되고, 부력항의 $\rho$ 를 절점값에서 쌍선형으로 보간하면 이쪽은 차수가 그대로다. 두 항이 서로 다른 함수공간에 있으므로 요소 내부에서 상쇄되지 않는다. 결과는 아무 일도 없어야 할 정수압 상태에서 요소마다 튀어나오는 가짜 속도이고, 밀도 대비가 클수록 커진다.

처방은 부력항의 밀도를 압력 기울기와 같은 차수로 떨어뜨려 평가하는 것이다. 개념적으로는 “상쇄해야 할 두 항을 같은 근사 차수에 두라”는 것이고, 같은 사상이 전산유체역학 쪽에서 well-balanced 도식으로 불린다. 정수압 균형을 이산 수준에서 정확히 만족시키지 못하는 도식은 정지 상태를 정지로 유지하지 못한다.

이 문제가 유한요소만의 것은 아니다. SEAWAT의 벤치마크 목록에 “회전 경계면” 문제가 들어 있는 이유가 정확히 이것 — 성층 밀도장을 회전시켜 놓고 가짜 속도가 나오는지 보는 시험이다. 밀도 결합 코드를 새로 짠다면 헨리 문제보다 먼저 돌려야 할 것이 정수압 성층 문제다. 답이 0이라 틀리면 바로 보인다.

5. 변포화 — SEAWAT에 없는 것[편집]

MODFLOW는 포화 유동 코드다. 자유면을 다루긴 하지만 그것은 “포화두께가 변한다”는 뜻이지 통기대의 모세관 물리를 푼다는 뜻이 아니고, UZF 패키지도 간이 저류 모형이지 리처즈 방정식이 아니다. SEAWAT은 이 제약을 그대로 물려받았다.

SUTRA는 $S_w(p)$ 와 $k_r(S_w)$ 를 넣어 불포화대를 직접 푼다. $S_w=1$ , $k_r=1$ 로 고정하면 위 유동식이 포화식으로 환원되므로, 포화 해석은 변포화 해석의 특수한 경우로 취급된다. 관심 물리가 다음 중 하나면 이 차이가 곧 코드 선택이 된다.

침투 지연과 모세관대. 지표 함양이 대수층에 도달하는 시간, 강우 이벤트에 대한 수위 반응.
해변 모세관대와 상부 염수 순환. 조석대 위쪽의 불포화 구간을 통과하는 염분 순환은 급격경계면 모형은 물론 포화 코드로도 못 본다.
토양 오염과 통기대 거동. 수송을 지표면부터 추적해야 하는 문제.
동토·계절 동결. 후기 판본의 동결-융해 확장(SUTRA-ICE 계열)이 얼음을 또 하나의 “불포화” 상태처럼 다뤄 영구동토 연구에 쓰인다.

대가도 분명하다. 불포화 유동은 $S_w(p)$ 와 $k_r$ 의 지독한 비선형 때문에 수렴이 어렵고, 건조한 초기 조건에서 습윤 전선이 들어올 때 시간스텝을 잘게 썰지 않으면 픽카드 반복이 발산한다. 그리고 반 게누흐텐 모형 같은 함수의 매개변수는 현장에서 거의 측정되지 않아, 물리를 더 정직하게 풀수록 측정되지 않은 매개변수가 늘어난다는 역설이 그대로 재현된다.

6. 결합과 시간 진행[편집]

한 시간스텝 안의 순서는 SEAWAT과 거의 같다. 현재 $U$ 로 $\rho$ · $\mu$ 를 계산하고 → 유동식을 풀어 $p$ 와 속도를 얻고 → 그 속도로 수송식을 풀어 새 $U$ 를 얻는다. 사용자가 반복 횟수를 지정하면 이 짝이 수렴할 때까지 반복되고(암시적 결합), 1이면 명시적 결합이다.

SUTRA에만 있는 손잡이가 하나 더 있다. 유동과 수송을 서로 다른 주기로 푼다. 압력해를 매 스텝이 아니라 $N$ 스텝마다 한 번만 갱신하도록 지정할 수 있는데(SUTRA 입력의 NPCYC·NUCYC), 배경 유동장이 거의 정상상태이고 농도만 천천히 움직이는 장기 모의에서 유동 솔버 비용을 통째로 아낀다. 유용하지만 위험한 옵션이기도 하다 — 밀도 결합이 강한 문제에서 이걸 켜면 부력 되먹임이 주기만큼 지연되어 결과가 조용히 틀린다. 켰다 껐을 때 답이 안 바뀌는지 확인하는 절차가 필수다.

7. 벤치마크[편집]

밀도 결합 코드의 검증 문화는 해수 침투 문서에 정리돼 있고, SUTRA는 그 문화를 상당 부분 만든 쪽이다.

문제	SUTRA가 여기서 보는 것
정수압 성층	정합 속도 — 가짜 속도가 0인가
헨리 문제	정상 염수 쐐기, 준해석해와 비교
엘더 문제	밀도 불안정 손가락, 격자 의존성
HYDROCOIN 암염돔	강한 밀도 대비의 순환
침투 기둥	불포화 습윤 전선(포화 코드에는 아예 없는 시험)

특히 엘더 문제를 지하수 분야의 표준 벤치마크로 만든 것이 SUTRA 계열의 계산이다. 원래 엘더(1967)의 실험은 다공층의 열대류였는데, 보스와 소자가 이것을 용질 버전으로 옮겨 코드 시험에 썼다. 그런데 이 문제는 격자를 조여도 하나의 답으로 수렴하지 않는다 — 손가락 개수와 중앙 상승/하강 여부가 격자와 솔버에 따라 갈리고, 이것은 코드 버그가 아니라 여러 안정해가 공존하는 분기 이론적 성질이다. 엘더 문제로 두 코드를 비교해 “어느 쪽이 맞다”고 결론 내리는 논문은 지금도 주기적으로 나오고 주기적으로 반박당한다.²

8. SEAWAT과의 대비[편집]

	SUTRA	SEAWAT
이산화	유한요소 + 질량집중	블록중심 유한차분
격자	사각형/육면체 요소망(변형 가능)	직교 구조격자
구조	단일 코드	MODFLOW + MT3DMS 결합
유동 미지수	압력	등가담수수두
불포화대	지원	미지원
수송 대상	용질 또는 열 1종	다중 화학종 + 열
반응	선형·프로인들리히·랭뮤어 흡착, 0·1차 생성소멸	MT3DMS의 반응 패키지 전체
생태계	얇다	MODFLOW 생태계 전부

요약하면 물리의 깊이는 SUTRA가, 주변부의 두께는 SEAWAT이 앞선다. 불포화 물리가 필요 없고 다중종 반응이 필요하며 PEST로 수천 번 돌릴 계획이라면 SEAWAT 또는 MODFLOW 6 쪽이 합리적이다. 반대로 통기대를 통과하는 염분 순환이나 지열 문제에서 요소망을 지층에 맞춰 깔아야 한다면 SUTRA 말고 대안은 FEFLOW 같은 상용 코드뿐이다.

생태계 이야기를 조금 더 하면, SUTRA의 오랜 약점은 격자를 만들기가 귀찮다는 것이었다. 유한요소 코드인데 메시 생성기가 딸려 오지 않아 한동안 Argus ONE 기반 SutraGUI와 SutraPlot에 의존했고, 나중에 ModelMuse가 SUTRA를 지원하면서 사정이 나아졌다. FloPy만큼 두터운 스크립트 층이 없다는 점은 여전히 대량 시나리오 작업에서 불리하게 작용한다.

9. 한계[편집]

화학이 얇다. 흡착 등온선과 1차 감쇠가 전부다. 다중종 반응망, 산화환원, 지구화학 평형은 범위 밖이며, 그쪽이 필요하면 PHREEQC 결합 코드나 OpenGeoSys 계열로 가야 한다.
한 번에 한 종. 용질과 열의 동시 수송(이중 확산 대류)은 본체가 못 한다. 파생판이 있지만 본류만큼 관리되지 않는다.
비선형 반복이 픽카드다. 밀도-농도 관계까지 야코비안에 넣는 완전 뉴턴-랩슨법이 아니라 고정점 반복이라, 강한 밀도 대비에서 수렴이 느리거나 진동한다. 완화계수를 손으로 만지는 작업이 남아 있다.
요소망이 축복이자 저주다. 지층에 맞춰 깔면 정확하지만, 심하게 뒤틀린 요소는 이류 지배 문제에서 정확도를 통째로 깎는다. 유한요소라고 격자 품질에서 자유로운 것은 전혀 아니다.
관측 자료 문제는 여전하다. 분산지수, 불포화 함수 매개변수, 비균질성 — 예측 오차를 지배하는 것은 이산화 오차가 아니라 이것들이라는 해수 침투의 결론이 SUTRA에서도 그대로다. 오히려 물리를 더 풀수록 미지의 매개변수가 늘어난다.³

그럼에도 SUTRA가 40년째 현역인 이유는 단순하다. 퍼블릭 도메인이고, 코드가 작고, 방정식이 문서에 전부 적혀 있다. 밀도 결합이나 변포화 알고리즘을 새로 시험해 보려는 사람이 참고 구현으로 집어 드는 코드가 아직 이것이고, 연구 논문의 파생판 개수로 치면 MODFLOW 계열 다음이다.

10. 관련 문서[편집]

해수 침투 · 지하수 유동 · 다공성 매질 유동 · 다르시 법칙 · 수리전도도
SEAWAT · MODFLOW · MT3DMS · MODPATH · FloPy · PEST
유한요소법 · 갤러킨 방법 · 혼합유한요소법 · 유한차분법
리처즈 방정식 · 대류-확산 방정식 · 페클레 수 · 수치분산
레일리-테일러 불안정 · 자연대류 · 분기 이론
뉴턴-랩슨법 · 희소행렬 · 검증 및 확인 · 불확실성 정량화

“통합 수송식” 아이디어는 우아하지만 딱 한 가지를 못 한다 — 이중 확산 대류(double-diffusive convection). 열과 염분이 확산속도가 달라서 생기는 층상 대류는 해양과 지열 분야의 실제 현상인데, 방정식이 하나뿐이면 원리적으로 표현이 안 된다. 하나의 식으로 둘을 푸는 설계는 “둘을 동시에는 못 푼다”는 제약과 동전의 양면이다. ↩
엘더 문제를 코드 검증용으로 쓰는 관행에는 근본적인 모순이 있다. 검증(verification)은 “알려진 답에 맞는가”를 보는 절차인데, 이 문제에는 알려진 답이 여러 개다. 그러니 엄밀히 말하면 엘더 문제는 검증 문제가 아니라 분기 구조 재현 시험이고, 두 코드가 다른 그림을 그렸다고 어느 한쪽이 틀렸다고 말할 수 없다. 물론 논문 심사에서 이 구분이 지켜지는 확률은 높지 않다. ↩
이 분야의 씁쓸한 국룰. 포화 모델을 돌리면 심의위원이 “불포화대는 왜 무시했느냐”고 묻고, 변포화 모델을 돌리면 “그 불포화 함수 매개변수는 어디서 났느냐”고 묻는다. 둘 다 정당한 질문이고 둘 다 답이 없다. 그래서 결국 민감도 해석 표를 두껍게 만들어 가는 것으로 합의가 이뤄진다. ↩