MODPATH

편집 역사 토론
소프트웨어 시뮬레이션 유체역학 마지막 수정: 2026-09-14 04:39:44

1. 개요[편집]

MODPATH
개발미국 지질조사국(USGS), 데이비드 폴록(D. W. Pollock)
초판1989
성격MODFLOW 유량장에서 입자 경로선을 뽑는 후처리기
핵심폴록의 준해석적 셀 통과 해
속도 보간각 성분을 자기 방향으로만 선형
물리이류만 — 분산·확산 없음
방향전진(오염원 → 하류) · 후진(관정 → 포획구역)
주요 판본3 · 5 · 6(2012) · 7(2016, MODFLOW 6 지원)

수두 등고선은 예쁘지만 아무도 그걸 보고 “이 오염물이 몇 년 뒤 우리 관정에 도착하는가”를 답하지 못한다. 그 질문에 답하려고 만든 것이 이 프로그램이다.

MODPATH는 MODFLOW가 계산한 셀 면 유량장을 읽어 가상 입자의 경로선(pathline)과 이동시간을 추적하는 USGS의 입자 추적 후처리기다. 1989년 데이비드 폴록이 만들었고, 그가 제안한 준해석적(semi-analytical) 셀 통과 방법이 지금도 이 프로그램의 전부라고 해도 된다. 관정 포획구역 설정, 오염 기원 추적, 대수층 체류시간 평가의 표준 도구.

MT3DMS와 자주 비교되는데 겨냥하는 질문이 다르다. MT3DMS는 농도를 계산하고 MODPATH는 물이 어디서 와서 어디로 가는가만 계산한다. 분산도 반응도 없고, 그래서 비용이 두세 자릿수 싸다. “농도까지는 필요 없고 물의 출처만 알고 싶다”는 요구가 생각보다 많고, 규제 문서에서 요구하는 것도 대개 그쪽이다.

2. 셀 유량에서 속도로[편집]

MODPATH가 읽는 것은 수두가 아니라 MODFLOW의 셀별 유량(cell-by-cell flow) 파일이다. 셀 면을 통과하는 유량 QQ 를 면적과 유효공극률로 나누면 그 면에서의 공극 유속이 나온다.

vx1=Qx1θ Δy Δz,vx2=Qx2θ Δy Δzv_{x1}=\frac{Q_{x1}}{\theta\,\Delta y\,\Delta z},\qquad v_{x2}=\frac{Q_{x2}}{\theta\,\Delta y\,\Delta z}

수두 기울기를 미분해서 속도를 만들지 않는다는 점이 중요하다. 수두를 미분하면 셀 면에서 값이 튀고 물수지와 어긋나는데, MODFLOW의 면 유량은 이미 셀 물수지를 정확히 만족하므로 그것을 그대로 쓰면 그 성질이 보존된다. MT3DMS가 같은 파일을 같은 이유로 읽는다.

여기서 공극률이 처음이자 마지막으로 등장한다. MODFLOW는 정상상태 유동에서 공극률을 전혀 쓰지 않는다 — 다르시 플럭스만 알면 되니까. 입자를 움직이려면 실제 물 알갱이의 속도가 필요하고, 그래서 MODPATH에서만 새로 요구되는 입력이 유효공극률이다. 그리고 이동시간은 여기에 정확히 반비례한다. 공극률을 0.3으로 놓느냐 0.15로 놓느냐가 10년 포획구역의 크기를 두 배로 바꾸는데, 이 값이 현장에서 측정되는 경우는 드물다. 포획구역 논쟁의 절반이 이 숫자에서 나온다.1

3. 폴록의 준해석적 해 — 이 문서의 킬포인트[편집]

입자 추적의 소박한 방법은 룽게-쿠타법으로 dx/dt=v(x)d\mathbf{x}/dt=\mathbf{v}(\mathbf{x}) 를 적분하는 것이다. 그러면 시간스텝을 골라야 하고, 스텝이 크면 셀을 건너뛰고, 작으면 느리다. 폴록은 그 선택 자체를 없앴다.

가정은 한 줄이다. 셀 안에서 각 속도 성분은 자기 좌표 방향으로만 선형으로 변한다.

vx(x)=vx1+Ax(x−x1),Ax=vx2−vx1Δxv_x(x)=v_{x1}+A_x\left(x-x_1\right),\qquad A_x=\frac{v_{x2}-v_{x1}}{\Delta x}

y,zy,z 도 똑같이. 그러면 xx 방향 운동방정식이 셋과 얽히지 않고 독립적으로 풀린다.

dxdt=vx1+Ax(x−x1)⟹vx(t)=vxp eAxt\frac{dx}{dt}=v_{x1}+A_x\left(x-x_1\right) \quad\Longrightarrow\quad v_x(t)=v_{xp}\,e^{A_x t}

즉 속도가 지수적으로 변하고, 위치는 닫힌 형태로 적분된다.

x(t)=x1+1Ax[vxp eAxt−vx1]x(t)=x_1+\frac{1}{A_x}\left[v_{xp}\,e^{A_x t}-v_{x1}\right]

vxpv_{xp} 는 현재 입자 위치의 속도다. 여기서 방향을 뒤집으면 면에 도달하는 시각을 해석적으로 얻는다. 입자가 x2x_2 면으로 나간다면

Δtx=1Axln⁡ ⁣vx2vxp\Delta t_{x}=\frac{1}{A_x}\ln\!\frac{v_{x2}}{v_{xp}}

여섯 개 면에 대해 이런 후보 시간을 전부 계산하고 가장 작은 것을 고르면, 그것이 실제 출구면과 통과시간이다. 그 시각을 위 위치식에 넣으면 출구점 좌표가 나온다.

결론: 입자는 셀 하나를 단 한 번의 계산으로 통과한다. 시간스텝이 없고, 셀 경계를 넘나드는 판정도 없고, 절단오차가 없다.2 Ax=0A_x=0 인 경우(그 방향 속도가 균일)는 지수식 대신 Δtx=(x2−xp)/vx\Delta t_x=(x_2-x_p)/v_x 로 처리하면 끝이다. 이것이 준해석적이라 불리는 이유이고, 40년 가까이 아무도 이 구성을 바꾸지 않은 이유다.

3.1. 왜 이 가정이 물수지와 정합적인가[편집]

“각 성분이 자기 방향으로만 선형”이라는 가정은 임의로 고른 것처럼 보이지만 그렇지 않다. 이 속도장의 발산을 계산해 보면

∇ ⁣⋅ ⁣v=Ax+Ay+Az=(Qx2−Qx1)+(Qy2−Qy1)+(Qz2−Qz1)θ Δx Δy Δz\nabla\!\cdot\!\mathbf{v}=A_x+A_y+A_z =\frac{\left(Q_{x2}-Q_{x1}\right)+\left(Q_{y2}-Q_{y1}\right)+\left(Q_{z2}-Q_{z1}\right)}{\theta\,\Delta x\,\Delta y\,\Delta z}

분자는 정확히 그 셀의 순 유출량이다. 원천·소멸이 없고 정상상태라면 MODFLOW의 물수지에 의해 0이므로, 보간된 속도장이 셀 안에서 정확히 발산 없는 장이 된다. 공짜로 얻어지는 성질이 아니다 — 삼선형 보간 같은 더 “정교한” 보간을 쓰면 이 성질이 깨지고, 입자가 밀폐된 영역에서 근거 없이 모이거나 흩어진다.

추가로, 면에 수직인 속도 성분은 이웃 셀과 같은 유량에서 나오므로 면을 가로질러 연속이다. 접선 성분은 불연속일 수 있지만 입자 궤적의 관통에는 영향이 없다. 즉 폴록의 구성은 정확도가 아니라 정합성을 위해 고른 것이고, 그래서 격자가 거칠어도 궤적이 물리적으로 말이 되는 방향으로 틀린다.

4. 전진·후진 추적과 포획구역[편집]

속도 부호를 뒤집으면 그대로 후진 추적이 된다. 이 단순한 사실이 MODPATH의 주 용도를 만든다.

  • 전진 추적 — 오염원 위치에 입자를 뿌리고 하류로 따라간다. “이 누출이 어느 관정에 도달하는가, 몇 년 뒤인가.”
  • 후진 추적 — 관정 셀에 입자를 채우고 거꾸로 따라간다. 입자들이 지표에 도달한 지점의 집합이 그 관정의 포획구역(capture zone)이다. 여기에 이동시간 등고선을 얹으면 2년·5년·10년 구역이 되고, 이것이 각국 상수원 보호구역 지정의 표준 근거다.

출력 형식도 이 용도에 맞춰 세 가지다. 종점(endpoint, 시작-끝만) · 경로선(pathline, 전체 궤적) · 시계열(time series, 지정 시각의 입자 위치 스냅숏). 포획구역 그림은 종점 파일에서, 애니메이션은 시계열에서 나온다.

비정상 유동에서는 경로선(pathline)·유적선(streakline)·유선(streamline)이 서로 다른 곡선이라는 유체역학 교과서의 구분이 여기서도 그대로 적용된다. MODPATH가 추적하는 것은 경로선이고, 정상 유동에서만 유선과 일치한다. 정상상태 모형의 유선망을 보고 비정상 문제의 이동시간을 논하는 것은 이 구분을 무시하는 것이다.

5. 약한 싱크 — 이 코드의 고전적 모호성[편집]

MODPATH 사용자가 반드시 한 번은 걸리는 문제. 어떤 셀에 우물이 있는데, 그 우물이 셀로 들어오는 물의 일부만 뽑아 간다고 하자. 나머지는 다음 셀로 흘러간다. 이 셀에 들어온 입자는 잡혀야 하는가, 통과해야 하는가?

물리적 답은 “확률적으로 일부만”이다. 그런데 입자는 쪼개지지 않는다. MODPATH는 이 모호성을 사용자에게 던진다.

  • 통과시킨다 — 포획구역이 과대평가된다(실제로 잡혔을 물까지 하류로 보내니 상류 추적 영역이 넓어진다).
  • 정지시킨다 — 포획구역이 과소평가된다. 약한 싱크에 걸린 입자가 거기서 끝나 버린다.
  • 임계값으로 나눈다 — 싱크가 셀 유입량의 몇 % 이상을 뽑을 때만 정지시킨다. 임계값은 또 하나의 손잡이다.

진짜 원인은 격자 해상도다. 셀이 우물 하나보다 훨씬 크기 때문에 “부분 포획”이라는 상태가 생긴다. 격자를 조이면 강한 싱크(들어온 물을 전부 뽑는 셀)로 수렴하고 모호성이 사라진다. 그래서 관정 주변 국소 세분화가 권장되며, 이것이 MODFLOW 6의 DISV 격자가 이 분야에서 환영받은 이유 중 하나다.

같은 종류의 모호성이 하천·배수·증발산 셀에서도 나온다. RIV 셀에 도달한 입자는 하천으로 나간 것인가 그 아래를 지나간 것인가? 답은 격자가 정하고, 보고서는 대개 이 선택을 적지 않는다.3

6. 무엇을 못 보나[편집]

MODPATH의 물리는 이류뿐이다. 이 한 줄에서 한계가 전부 나온다.

  • 분산·확산이 없다. 포획구역 경계가 날카로운 선으로 그려지지만, 실제 오염운의 가장자리는 퍼져 있다. MODPATH의 10년 선 바로 바깥이 안전하다는 뜻이 결코 아니다. 확률적 경계를 원하면 무작위 변위를 더한 무작위 행보 입자추적(RWPT)이나 MT3DMS로 가야 한다.
  • 반응이 없다. 흡착과 붕괴는 표현되지 않는다. MODPATH는 구역별 지연계수를 받아 이동시간을 그 배수로 늘려 주는 정도의 기능만 제공하며, 이것은 균일 선형 흡착의 아주 거친 근사다.
  • 밀도가 균일하다고 본다. 염수 쐐기 위에서 입자를 추적하면 부력을 무시한 궤적이 나온다. 해수 침투 문제에서 MODPATH 결과를 그대로 쓰면 안 되는 이유.
  • 유동해가 틀리면 경로선도 틀린다. 그리고 이쪽 오차가 압도적으로 크다. 이동시간은 공극률과 수리전도도에 직접 비례·반비례하므로, 이 두 값의 불확실성이 그대로 곱해져서 나온다. 포획구역을 하나의 선으로 보고하는 관행 자체가 불확실성 정량화 관점에서는 이미 문제이며, 요즘은 PEST 계열로 만든 매개변수 앙상블마다 MODPATH를 돌려 포획 확률 지도를 그리는 쪽이 권장된다.

7. 계보[편집]

  • MODPATH 1~3. 폴록의 원 구성이 확립된 시기. 3에서 MODFLOW-96 파일 형식과 비정상 유동 추적이 정리됐다.
  • MODPATH 5. 오래 쓰인 안정 판본. 상용 GUI 다수가 이 형식에 묶여 있다.
  • MODPATH 6 (2012). 입력 구조를 정리하고 추적 옵션(약한 싱크 처리, 정지 구역)을 명시적으로 노출했다.
  • MODPATH 7 (2016). MODFLOW 6과 MODFLOW-USG 지원이 핵심. 구조격자(DIS)뿐 아니라 평면 다각형 격자(DISV)까지 다룬다.

7의 비구조 격자 지원은 생각보다 큰 일이었다. 폴록의 구성은 직육면체 셀 + 마주 보는 면 쌍이라는 기하에 기대고 있어서, 육각형이나 보로노이 다각형 셀에서는 ”xx 방향 두 면” 자체가 정의되지 않는다. 다각형 셀을 부분 영역으로 쪼개 각 조각에서 같은 구성을 되살리는 식의 처리가 필요하고, 그 과정에서 원래의 발산 정합성을 어떻게 유지하는가가 설계의 핵심 쟁점이 된다. “격자를 유연하게 만들면 후처리기가 고생한다”는 법칙이 여기서도 성립한다.

곁가지로 MODPATH-OBS가 있다. 추적 결과에서 이동시간이나 농도 비슷한 양을 관측값으로 뽑아 주는 도구인데, 이것이 있으면 입자 추적 결과 자체를 보정 자료로 쓸 수 있다. 삼중수소 연대나 관정 혼합 연령 같은 자료를 PEST에 먹일 수 있게 되는 것이라, 수두만으로는 제약되지 않는 공극률·연직 흐름 분배를 잡는 데 쓰인다.

8. 관련 문서[편집]

9. Footnotes[편집]

  1. 더 나쁜 것은 이 값이 보정으로도 잘 잡히지 않는다는 점이다. 정상상태 수두 자료에는 공극률이 아예 들어가지 않으므로, 수두만 맞춰서 보정한 모형은 공극률에 대해 아무 정보도 갖고 있지 않다. 그런데도 그 모형으로 10년 포획구역을 그려 제출한다. 공극률을 제약하려면 연대 추적자(삼중수소, CFC, 헬륨-4)나 양수시험 추적자 자료가 필요하고, 그런 자료가 있는 현장은 드물다. ↩

  2. “절단오차가 없다”는 말을 “정확하다”로 읽으면 곤란하다. 가정된 속도장 안에서 궤적을 정확히 푼 것일 뿐, 그 속도장 자체가 실제 속도장의 거친 근사다. 즉 오차는 시간적분이 아니라 공간 보간에 통째로 들어 있다. 격자를 조이면 궤적이 좋아지는 것은 그래서이며, 룽게-쿠타법으로 스텝을 아무리 잘게 썰어도 같은 속도장을 쓰는 한 답은 거의 같은 곳으로 간다. 시간적분의 정확도를 올리는 데 시간을 쓰는 것은 이 문제에서 헛수고다. ↩

  3. 약한 싱크 처리 옵션을 바꿔 가며 세 번 돌려 보면 포획구역 면적이 몇십 퍼센트씩 달라지는 일이 흔하다. 정직한 보고서라면 그 세 개를 다 싣고 “이 폭이 격자 해상도에서 온 것”이라고 적어야 한다. 실제로는 셋 중 제일 그럴듯한 그림 하나가 실린다. 어느 것이 제일 그럴듯한가는 대개 그 그림이 누구의 부지를 포함하느냐로 결정된다. ↩