MT3DMS

1. 개요[편집]

MT3DMS Modular 3-D Multi-Species Transport Model
개발	정춘묘(Chunmiao Zheng) 외 — USEPA · 미 육군 공병단 지원
계보	MT3D(1990) → MT3DMS(1999) → MT3D-USGS(2016)
역할	MODFLOW 유동해를 받아 용질 수송만 푸는 후속 코드
미지수	셀 중심 농도 C (화학종별)
이류 해법	FD · MOC · MMOC · HMOC · TVD(ULTIMATE)
솔버	GCG(일반화 켤레기울기) — 분산·반응·원천항 음해
제약	격자 페클레 수 · 쿠랑 수
언어·라이선스	Fortran · 퍼블릭 도메인

물이 어디로 흐르는지는 MODFLOW가 답한다. 그 물에 섞인 것이 어디로 가는지는 다른 코드가 답한다. 그 다른 코드가 20년째 이것이다.

MT3DMS는 MODFLOW가 계산한 셀 면 유량장을 입력으로 받아 다공성 매질에서의 이류·분산·화학반응 수송을 푸는 모듈형 용질 수송 코드다. 이름의 MS는 다중 화학종(multi-species)을 뜻하며, 1990년 정춘묘가 만든 MT3D를 1998~99년에 대폭 확장한 것이 현재 계보의 뿌리다. 흔한 오해와 달리 USGS 코드가 아니다 — 미 환경보호청(USEPA)과 육군 공병단의 지원으로 개발됐고, USGS가 관리하는 것은 뒤에 나온 파생판 MT3D-USGS다.¹

푸는 방정식은 대류-확산 방정식에 반응항을 붙인 것뿐이라 새로울 게 없다. 이 문서의 주인공은 방정식이 아니라 이류항을 어떻게 이산화할 것인가라는 선택지다. MT3DMS는 그 선택지를 다섯 개나 사용자에게 열어 두었고, 그 다섯 개가 곧 “수치분산이냐 인위진동이냐”라는 수치 소산과 분산의 고전적 딜레마를 그대로 메뉴로 만든 것이다.

2. 무엇을 푸는가[편집]

화학종 $k$ 에 대해

$$\frac{\partial\left(\theta C^k\right)}{\partial t} =\nabla\!\cdot\!\left(\theta\,\mathbf{D}\cdot\nabla C^k\right) -\nabla\!\cdot\!\left(\theta\,\mathbf{v}\,C^k\right) +q_s C_s^k+\sum_n R_n$$

$\theta$ 는 공극률, $\mathbf{v}$ 는 공극 유속(다르시 플럭스를 $\theta$ 로 나눈 것), $q_s$ 는 단위 체적당 원천·소멸 유량, $\mathbf{D}$ 는 분산 텐서다. 분산 텐서는 셰이데거-베어 형태로 종분산지수 $\alpha_L$ 과 횡분산지수 $\alpha_{TH},\alpha_{TV}$, 그리고 분자확산 $D^*$ 로 구성된다.

반응항 쪽에서 실무적으로 제일 많이 쓰이는 것은 선형 평형 흡착이고, 이때 좌변에 지연계수가 붙는다.

$$R=1+\frac{\rho_b}{\theta}\frac{\partial\bar{C}}{\partial C} \quad\xrightarrow{\ \bar{C}=K_dC\ }\quad R=1+\frac{\rho_b K_d}{\theta}$$

효과는 단순하다 — 전선이 물보다 $R$ 배 느리게 간다. 프로인들리히·랭뮤어 비선형 흡착 등온선, 1차 감쇠(모종 사슬 붕괴 포함), 그리고 이중영역 물질전달(dual-domain mass transfer)까지 RCT 패키지에 들어 있다. 마지막 것이 MT3D 대비 MT3DMS의 중요한 추가인데, 이동영역과 비이동영역 사이의 1차 교환을 두어 관측 파과곡선의 긴 꼬리를 재현한다. 비균질 매질에서 꼬리가 안 맞을 때 분산지수를 키워 억지로 맞추는 대신 쓸 수 있는 물리적으로 더 정직한 손잡이다.

패키지 구조는 MODFLOW를 그대로 본떴다. BTN(기본 수송) · ADV(이류) · DSP(분산) · SSM(원천·소멸 혼합) · RCT(반응) · GCG(솔버). 파일 하나 = 물리 하나라는 조립 방식도 같다.

3. 유동해를 어떻게 받는가[편집]

MT3DMS는 유동을 풀지 않는다. MODFLOW 쪽에서 LMT(Link-MT3DMS) 패키지를 켜면 셀별 수두·포화두께·셀 면 유량·원천항이 담긴 FTL 연결 파일이 떨어지고, MT3DMS는 그것만 읽는다. 이 설계가 만드는 성질이 몇 가지 있다.

• 격자가 동일해야 한다. 유동 격자와 수송 격자가 같은 행·열·층이다. 유동은 거친 격자로 충분한데 수송은 전선 근처에서 조밀해야 하는 상황이 흔한데, 이 구조에서는 유동 쪽을 함께 조여야 한다. 수송 문제의 계산비용이 유동보다 한두 자릿수 비싼 주된 이유.
• 질량 보존이 공짜로 따라온다. MODFLOW의 면 유량은 이미 셀 물수지를 정확히 만족하므로, 그 유량으로 만든 속도장은 이산 수준에서 발산이 0이다. 수송 코드가 속도장을 스스로 보간해서 만들었다면 이 성질이 깨진다 — MODPATH가 같은 자료를 같은 이유로 읽는다.
• 단방향 결합이다. 농도가 유동에 영향을 주지 않는다는 가정이 들어 있다. 밀도가 농도에 의존하면 이 가정이 깨지고, 그때 필요한 것이 유동-수송 되먹임 루프를 건 SEAWAT이다.

시간 축도 이중 구조다. MODFLOW의 유동 시간스텝 하나 안에서 MT3DMS는 자신의 안정 조건에 따라 수송 스텝을 여러 개로 쪼갠다. 유동 스텝이 30일이어도 수송 스텝은 며칠이 되는 일이 흔하고, 그래서 총 실행시간의 대부분이 수송 쪽에서 소모된다.

4. 이류항 — 이 문서의 본론[편집]

농도 전선은 본질적으로 불연속에 가깝다. 그 불연속을 격자 위에 올릴 때 무엇을 포기할 것인가가 전부다.

4.1. 1. 표준 유한차분 (MIXELM = 0)[편집]

상류가중이면 무조건 안정하고 질량이 보존되며 농도가 음수로 안 간다. 대가는 **수치분산**이다. 1차 상류차분의 절단오차는 $\frac{1}{2}v\Delta x$ 크기의 인위적 분산계수와 같으므로, 이것을 실제 분산 $v\alpha_L$ 과 비교하면 조건이 바로 나온다.

$$\frac{\Delta x}{2}\lesssim\alpha_L \quad\Longleftrightarrow\quad \mathrm{Pe}=\frac{\Delta x}{\alpha_L}\lesssim 2$$

격자 페클레 수가 2를 넘으면 계산된 전이대 두께를 물성이 아니라 격자가 정한다. 중앙가중으로 바꾸면 수치분산은 사라지지만 같은 $\mathrm{Pe}>2$ 에서 진동이 시작되고, 진동이 농도를 음수로 만든다. 음의 농도는 그 자체로도 곤란하지만 SEAWAT처럼 밀도 결합이 걸려 있으면 없던 부력을 만들어내는 진짜 재앙이 된다.

현실에서 $\alpha_L$ 은 미터 단위인데 모형 도메인은 킬로미터라, $\mathrm{Pe}\le2$ 를 지키려면 격자가 감당 못 할 크기가 된다. MT3DMS 매뉴얼조차 실무에서는 4 정도까지 눈감아 주는 편이라고 적어 둘 정도다. 이 구조적 곤란이 나머지 네 개 옵션이 존재하는 이유다.

4.2. 2~4. 특성선 계열 — MOC · MMOC · HMOC[편집]

MOC(method of characteristics)는 이류를 격자가 아니라 입자로 옮긴다. 셀마다 입자를 여럿 뿌려 각자 농도를 지니게 하고, 속도장을 따라 이동시킨 뒤(전진 추적), 셀에 들어온 입자들의 평균으로 농도를 갱신한다. 이류만 담당하고 분산·반응은 유한차분이 이어받는다. 격자 확산이 원리적으로 0이라 급격한 전선을 날카롭게 유지하는 것이 최대 장점.

단점이 셋이다. 입자를 많이 뿌려야 해서 메모리가 무겁고, 셀별 입자 수가 들쭉날쭉해 농도장이 얼룩덜룩(numerical lumping)해지며, 결정적으로 질량 보존이 정확하지 않다. 입자 삽입·삭제와 평균화 과정에서 질량이 새는데, 물수지표의 불일치율이 몇 퍼센트씩 나오는 것을 보고 처음 당황하는 지점이다.

MMOC(modified MOC)는 반대로 간다. 셀 중심마다 입자를 하나씩 놓고 시간스텝만큼 역방향으로 추적해 상류 위치를 찾은 뒤, 그 위치의 농도를 보간해 가져온다. 입자 수가 셀 수로 고정되니 싸고 얼룩도 없다. 대신 상류 위치에서 농도를 보간하는 순간 다시 수치분산이 들어온다 — 전선이 급격할수록 심하다.

HMOC(hybrid)는 둘을 자동으로 섞는다. 농도 기울기가 임계값을 넘는 셀 주변에서만 MOC를 쓰고 나머지는 MMOC로 처리한다. 임계값은 사용자가 준다. 개념적으로는 적응형 기법의 정석이고, 실제로는 임계값이 또 하나의 손잡이로 남는다.

4.3. 5. TVD — ULTIMATE (MIXELM = -1)[편집]

레너드의 ULTIMATE 제한자를 3차 상류 도식에 얹은 유한체적 정식화다. 전선 근처에서 국소적으로 차수를 떨어뜨려 새로운 극값이 생기지 않도록(전변분 감소) 플럭스를 제한한다. MUSCL이나 전산유체역학 쪽 고해상도 도식과 정확히 같은 계보이며, 진동 없이 질량을 보존하고 수치분산도 1차 상류차분보다 훨씬 작다.

값은 있다. 명시적 도식이라 쿠랑 수 제약을 받는다.

$$\mathrm{Cr}=\frac{|v|\,\Delta t}{\Delta x}\le 1$$

우물 근처처럼 속도가 큰 셀 하나가 전체 수송 스텝을 지배해 스텝 수가 폭증하는 상황이 흔하다. 그래도 오늘날의 기본 권고는 일단 TVD다. MOC 계열의 질량 누수와 손잡이를 감수할 이유가 대부분의 문제에서 없기 때문.

4.4. 정리[편집]

옵션	수치분산	진동	질량보존	시간 제약	비용
FD 상류	큼	없음	정확	음해 가능	싸다
FD 중앙	작음	Pe > 2 에서 발생	정확	음해 가능	싸다
MOC	거의 없음	얼룩	부정확	Cr 제약	비쌈
MMOC	중간	없음	양호	Cr 제약	중간
HMOC	작음	없음	양호	Cr 제약	중간
TVD	작음	없음	정확	Cr ≤ 1	중간

5. GCG 솔버 — 무엇이 음해로 풀리나[편집]

MT3DMS의 또 하나의 기여가 GCG(Generalized Conjugate Gradient) 솔버다. 분산·원천소멸·반응 항을 전부 음해로 묶어 켤레기울기법 계열로 푼다. 이게 없던 MT3D 시절에는 분산항도 명시적이라 다음 조건에 묶여 있었다.

$$\Delta t\le\frac{0.5}{\dfrac{D_{xx}}{\Delta x^2}+\dfrac{D_{yy}}{\Delta y^2}+\dfrac{D_{zz}}{\Delta z^2}}$$

격자를 반으로 줄이면 허용 스텝이 4분의 1이 되는 전형적인 확산 제약이다. GCG가 이 제약을 통째로 걷어내면서 수송 모델링의 격자 해상도 한계가 크게 완화됐다.

주의할 것은 이류는 여전히 별도라는 점이다. FD 이류를 음해로 돌리면 스텝 제약이 없어지지만 수치분산을 감수해야 하고, MOC·TVD를 쓰면 그쪽은 명시적이라 쿠랑 제약이 살아 있다. “GCG를 켰는데 왜 스텝이 안 커지느냐”는 질문의 답은 거의 항상 여기 있다.

6. 계보[편집]

• MT3D (1990). 정춘묘가 USEPA 지원으로 개발. MOC 계열 3종이 이때 이미 들어 있었다.
• MT3DMS (1998~1999). 다중 화학종 구조, TVD(ULTIMATE), GCG 음해 솔버, 이중영역 물질전달. 육군 공병단 보고서로 문서화됐고 사실상 지금 쓰이는 판본이다.
• SEAWAT (2002~2008). MT3DMS를 MODFLOW와 밀도 결합 루프로 묶어 한 실행 파일에 넣은 것. 열 수송은 온도를 화학종 하나로 취급해 구현했다.
• MT3D-USGS (2016). USGS가 이어받아 재작성. MODFLOW-NWT와 물리기반으로 연결되고, 하천(SFT)·호소(LKT)·불포화대(UZT) 수송, 오염 처리 시스템(CTS) 모의가 추가됐다. MODFLOW가 UZF·SFR로 표현한 것을 수송 쪽에서도 일관되게 받는 것이 핵심 동기.
• MODFLOW 6 GWT. MODFLOW 6 본체가 지하수 수송 모형을 들이면서, 예전에 외주 주던 일이 같은 실행 파일 안으로 들어왔다. 비구조 격자와 뉴턴 정식화가 딸려 오고, 유동-수송 결합이 예외가 아니라 일반 교환(exchange) 기능의 한 사례가 된다.

그럼에도 MT3DMS가 아직 현역인 이유는 SEAWAT과 같다 — 20년치 모형 자산과 인허가 문서가 이 형식으로 존재하고, FloPy가 계속 지원하며, 논문의 벤치마크가 이 결과를 기준으로 인용된다.

7. 실무에서의 현실[편집]

• 수송 보정은 유동 보정보다 훨씬 나쁜 역문제다. 농도 관측은 수두 관측보다 드물고, 검출한계 이하 자료가 태반이며, 시료 채취 시점이 제각각이다. 그래서 PEST로 $\alpha_L$ 과 $K_d$ 를 동시에 추정하면 둘이 서로를 대신 흡수해 버린다.
• 분산지수는 물성이 아니라 손잡이다. 척도 의존성이 심해 실험실 값과 현장 값이 두 자릿수 차이 난다. 전이대가 안 맞으면 제일 먼저 손대는 값이고, 그래서 이 값에 물리적 의미를 부여하는 순간 논쟁이 시작된다.
• 격자를 조이는 것이 항상 답은 아니다. 페클레 조건을 만족시키려고 격자를 조이면 쿠랑 조건이 스텝을 조이고, 실행 시간이 세제곱 이상으로 늘어난다. 예산 안에서는 대체로 TVD + 적당한 격자 + 민감도 해석이 현실적 타협이다.
• 물수지를 두 개 다 본다. MODFLOW의 물 물수지와 MT3DMS의 용질 물수지는 별개다. 물 쪽이 0.1%인데 용질 쪽이 5%인 상황이 MOC에서는 예사로 나오고, 그 표를 안 보고 넘어간 결과가 보고서에 그대로 실린다.²
• 그리고 이 바닥의 오래된 농담. “전선이 너무 퍼졌는데요?” — 그러면 분산지수를 줄이거나, 격자를 조이거나, 둘 다 안 하고 그래프의 세로축 범위를 바꾼다.³

8. 관련 문서[편집]

• MODFLOW · SEAWAT · SUTRA · MODPATH · FloPy · PEST
• 대류-확산 방정식 · 페클레 수 · CFL 조건
• 수치분산 · 수치 소산과 분산 · 전변분 감소 · MUSCL · 특성곡선법
• 지하수 유동 · 해수 침투 · 다공성 매질 유동 · 수리전도도
• 유한차분법 · 유한체적법 · 켤레기울기법 · 희소행렬
• 역문제 · 민감도 해석 · 불확실성 정량화 · 검증 및 확인

9. Footnotes[편집]

1. 이 혼동이 널리 퍼진 데는 이유가 있다. MT3DMS가 MODFLOW와 한 몸처럼 쓰이고, USGS가 배포하는 SEAWAT 안에 통째로 들어가 있으며, 결국 MT3D-USGS가 나왔기 때문이다. 그래도 원 개발자와 자금원을 정확히 적는 것은 라이선스 문제가 아니라 예의의 문제다 — 퍼블릭 도메인이라 아무나 가져다 팔 수 있는 것과, 누가 만들었는지를 적는 것은 별개다. ↩

2. 용질 물수지가 크게 어긋났을 때 흔한 응급처치가 입자 수를 늘리는 것인데, 이것은 증상을 완화할 뿐이다. MOC의 질량 누수는 입자 삽입·삭제 규칙에서 오는 구조적 성질이라 입자를 늘려도 0으로 가지 않는다. 진짜 처방은 TVD로 옮기는 것이고, 그러면 대개 실행시간이 늘어난다. 정확한 답과 빠른 답 사이의 선택이 또 나온다. ↩

3. 물론 진짜로 해야 할 일은 셋 중 어느 것도 아니다. 계산한 전이대 두께가 $\alpha_L$ 이 만든 것인지 $\Delta x$ 가 만든 것인지 구분하는 것이다. 방법은 단순하다 — 분산지수를 0으로 놓고 한 번 돌려 본다. 그래도 전선이 퍼져 있으면 그 퍼짐은 전부 격자가 만든 것이다. 30분이면 끝나는 이 시험을 안 하고 보고서를 쓰는 경우가 놀랄 만큼 많다. ↩