리처즈 방정식

1. 개요[편집]

리처즈 방정식 Richards equation
발표	L. F. Richardson (1922) · L. A. Richards (1931)
대상	불포화대(vadose zone) 물 이동
혼합형	∂θ/∂t = ∇·[K(h)(∇h + ez)]
보유곡선	반 헤누흐텐-무알렘 · 브룩스-코리
성격	퇴화 비선형 포물형 (포화 시 타원형)
표준 이산화	Celia(1990) 혼합형 + 수정 피카르
대표 코드	HYDRUS · ParFlow · SWAP

미지수가 하나뿐인 스칼라 방정식인데, 수치해석 학회에서 이 방정식 얘기가 나오면 다들 눈을 피한다.

리처즈 방정식(Richards equation)은 공극이 물과 공기로 함께 채워진 불포화 다공성 매질에서, 기체상의 압력을 대기압으로 고정한 채 물의 이동만을 기술하는 지배 방정식이다. 혼합형(mixed form)으로 쓰면

$$\frac{\partial\theta}{\partial t}=\nabla\cdot\!\left[K(h)\left(\nabla h+\mathbf{e}_z\right)\right]+S$$

이고, $\theta$ 는 체적함수량, $h$ 는 압력수두(불포화대에서 음수), $K(h)$ 는 불포화 수리전도도, $\mathbf{e}_z$ 는 연직 상향 단위벡터(중력항), $S$ 는 뿌리 흡수 같은 싱크항이다.

물리적 내용물은 단순하다. 불포화 상태의 다르시 법칙에 질량보존을 붙인 것이 전부다. 포화 지하수 방정식과 달라지는 지점은 딱 하나 — 공극에 공기가 있으면 물이 지나갈 통로가 줄고 모관력이 물을 붙잡으므로, $K$ 와 $\theta$ 가 $h$ 의 강한 비선형 함수가 된다. 이 한 가지가 방정식 전체의 성격을 바꾸고, 수치해석 쪽 악명 전부의 원인이 된다.¹

강우 침투, 관개 설계, 농약·질산염 지하수 오염, 사면 안정, 매립장 차수층, 지표수-지하수 결합 모델의 연결 고리가 전부 이 식이다. 계산이 감당 안 될 때 쓰는 축약형이 그린-암프트 모형이고, 그쪽은 이 방정식의 급격전선 극한이다.

2. 세 가지 정식화[편집]

같은 물리인데 미지수를 무엇으로 잡느냐에 따라 세 얼굴이 나온다. 이게 그냥 취향 문제가 아니라 질량보존이 지켜지느냐 마느냐를 가른다.

압력수두형(h-form). 연쇄법칙으로 $\partial\theta/\partial t=C(h)\,\partial h/\partial t$ 를 쓴다. $C(h)=\mathrm{d}\theta/\mathrm{d}h$ 가 비수분용량(specific moisture capacity)이다.

$$C(h)\frac{\partial h}{\partial t}=\nabla\cdot\!\left[K(h)(\nabla h+\mathbf{e}_z)\right]$$

장점: 포화·불포화 영역을 한 미지수로 다룰 수 있고, 층상 토양의 물질 경계에서 $h$ 는 연속이라 경계조건이 자연스럽다. 단점: 뒤에서 볼 질량 손실.

함수량형(θ-form). $D(\theta)=K\,\mathrm{d}h/\mathrm{d}\theta$ (토양수분확산계수)를 도입하면 이류-확산 꼴이 된다.

$$\frac{\partial\theta}{\partial t}=\nabla\cdot\!\left[D(\theta)\nabla\theta\right]+\frac{\partial K}{\partial z}$$

대류-확산 방정식과 같은 골격이라 수치적으로 다루기 편하고 질량보존도 자동이다. 하지만 쓸 수 없는 경우가 너무 많다. 포화되면 $\mathrm{d}h/\mathrm{d}\theta\to\infty$ 라 $D$ 가 발산하고, 층상 토양에서는 물질 경계에서 $\theta$ 가 불연속이라(모래와 점토가 같은 $h$ 에서 전혀 다른 $\theta$ 를 갖는다) $\nabla\theta$ 가 의미를 잃는다. 균질 불포화 1차원 해석용 형태로 봐야 한다.

혼합형(mixed form). 시간항은 $\theta$ 로, 공간항은 $h$ 로 두고 미지수는 $h$ 로 푼다. 위 개요의 식이 이것이다. 두 형태의 장점만 가져온 대신, $\theta$ 와 $h$ 를 잇는 이산화를 제대로 해야 한다는 숙제가 남는다. 오늘날 표준이다.

3. 보유곡선 — 비선형의 정체[편집]

$\theta(h)$ 를 토양수분보유곡선(water retention curve), $K(\theta)$ 또는 $K(h)$ 를 불포화 전도도 함수라 한다. 이 두 개가 리처즈 방정식의 물성 입력 전부다.

반 헤누흐텐(1980) 모형이 압도적 표준이다. 유효포화도 $S_e=(\theta-\theta_r)/(\theta_s-\theta_r)$ 에 대해

$$S_e=\left[1+(\alpha\lvert h\rvert)^{n}\right]^{-m},\qquad m=1-\frac1n$$

이고, 여기에 무알렘(1976) 의 공극연결 모형을 결합하면 전도도가 닫힌 형태로 나온다.

$$K(S_e)=K_s\,S_e^{\,\ell}\left[1-\left(1-S_e^{1/m}\right)^{m}\right]^{2},\qquad \ell=0.5$$

$m=1-1/n$ 이라는 구속을 건 이유가 이것 — 그래야 무알렘 적분이 닫힌 형태로 떨어진다. 이론적 필연이 아니라 적분이 풀리게 하려고 건 제약이다.²

브룩스-코리(1964) 는 더 오래된 멱함수형이다. 공기침입수두 $h_b$ 와 공극크기분포지수 $\lambda$ 로

$$S_e=\left(\frac{h_b}{\lvert h\rvert}\right)^{\lambda}\ (\lvert h\rvert>\lvert h_b\rvert), \qquad K=K_s\,S_e^{(2+3\lambda)/\lambda}$$

$h_b$ 에서 기울기가 꺾이는 비매끄러운 곡선이라 뉴턴 반복이 싫어하지만, 물리적으로는 “공기침입압 이하에서는 물이 안 빠진다”는 사실을 정직하게 표현한다. 반 헤누흐텐은 그 꺾임을 매끄럽게 뭉갠 대신 미분가능성을 얻었다.

카슬-패리시(1988)가 정리해 코드 기본값으로 굳은 대표값을 보면 비선형의 규모가 보인다.

토양	θs	α (1/cm)	n	Ks (cm/d)
모래	0.43	0.145	2.68	713
양토	0.43	0.036	1.56	25
미사	0.46	0.016	1.37	6
점토	0.38	0.008	1.09	4.8

포화 전도도가 150배 차이나는 것은 놀랍지도 않다. 진짜 문제는 하나의 토양 안에서 $K$ 가 움직이는 폭이다. 모래에서 $S_e$ 가 1에서 0.3으로 떨어지면 $K$ 는 대여섯 자릿수 떨어진다. 습윤전선 하나를 사이에 두고 셀과 이웃 셀의 계수가 $10^5$ 배 차이나는 상황이 매 스텝 벌어진다는 뜻이며, 이게 이 방정식이 유명한 진짜 이유다.³

4. 왜 그렇게 안 풀리는가[편집]

리처즈 방정식은 퇴화(degenerate) 비선형 포물형 방정식으로 분류된다. “퇴화”가 붙은 이유가 곧 고통의 목록이다.

• 포화되면 포물형이 아니게 된다. $h\ge0$ 에서 $C(h)=\mathrm{d}\theta/\mathrm{d}h=0$ 이라 시간 미분항이 사라지고 방정식이 타원형으로 바뀐다. 행렬의 대각 우세가 무너지고, 시간 전진이 아니라 제약조건을 푸는 문제가 된다. 실무 코드는 포화 영역에 작은 비저류계수 $S_s$ 를 넣어 퇴화를 인위적으로 없앤다.
• 아주 마르면 반대쪽으로 퇴화한다. $K\to0$ 이면 확산이 사라져 전파가 멈춘다. 초기조건이 극도로 건조한 문제(사막 토양에 물을 붓는 고전 벤치마크)에서 반복이 죽는 이유.
• 습윤전선은 거의 불연속이다. 전선 두께가 수 mm~cm 수준이라, 격자가 이를 해상하지 못하면 계수 평균이 통째로 틀린다. 셀 경계 $K$ 를 산술평균으로 잡으면 전선이 실제보다 빨리 가고, 조화평균으로 잡으면 느리게 간다. 기하평균이나 상류가중(upstream weighting)을 쓰는 이유이며, 이 선택 하나로 침투 깊이가 수십 % 달라진다.
• 비단조 응답. $C(h)$ 는 종 모양이라 전선 부근에서 몇 자릿수로 치솟았다 꺼진다. 야코비안이 스텝마다 딴판이 된다는 뜻.

이 모든 것이 합쳐진 결과가 시간 간격 붕괴다. 습윤전선이 격자를 하나 통과할 때마다 비선형 반복이 안 붙고, 코드는 $\Delta t$ 를 반으로 자르고 다시 시도한다. 밤새 돌린 연속 모의가 “10년 중 3개월”에서 멈춰 있는 광경은 이 바닥의 국룰이다.

5. Celia 1990 — 질량이 새는 이유[편집]

1990년 이전 코드는 대부분 h-form을 후진 오일러로 이산화했다. 반복 $m$ 단계에서

$$C(h^{n+1,m})\,\frac{h^{n+1,m+1}-h^{n}}{\Delta t}=\nabla\cdot\!\left[K^{m}(\nabla h^{n+1,m+1}+\mathbf{e}_z)\right]$$

이 형태의 치명적 결함은 저장 변화량을 $C\,\Delta h$ 로 근사한다는 것이다. 실제 저장 변화는 $\theta(h^{n+1})-\theta(h^n)$ 인데, $C$ 가 $h$ 에 따라 급격히 변하므로 이 근사의 오차는 반복을 아무리 돌려도 사라지지 않는다. 반복은 자기 자신이 정의한 잘못된 방정식에 수렴한다. 결과가 그 유명한 질량 손실 — 침투 문제에서 총 유입량 대비 몇 %에서 수십 %까지 물이 증발하듯 사라지고, $\Delta t$ 를 줄여야만 줄어든다.

Celia, Bouloutas, Zarba(1990)의 처방은 혼합형을 그대로 이산화하고, $\theta$ 만 테일러 전개하는 것이다.

$$\theta^{n+1,m+1}\approx\theta^{n+1,m}+C^{n+1,m}\left(h^{n+1,m+1}-h^{n+1,m}\right)$$

이걸 혼합형에 넣으면 미지수를 증분 $\delta=h^{n+1,m+1}-h^{n+1,m}$ 로 하는 선형계가 나온다.

$$\frac{C^{m}}{\Delta t}\delta-\nabla\cdot\!\left[K^{m}\nabla\delta\right] =\nabla\cdot\!\left[K^{m}(\nabla h^{m}+\mathbf{e}_z)\right]-\frac{\theta^{m}-\theta^{n}}{\Delta t}$$

핵심은 우변이다. 저장항이 $C\Delta h$ 가 아니라 $\theta^m-\theta^n$ 이라는 참값으로 들어가 있다. 반복이 수렴해 $\delta\to0$ 이 되면 남는 것은 정확한 질량보존식이므로, 질량균형은 $\Delta t$ 의 크기와 무관하게 반복 수렴 오차 수준으로 지켜진다. 이 한 줄짜리 재배치가 이후 30년의 표준이 됐고, 오늘날 HYDRUS를 포함한 사실상 모든 코드가 이 수정 피카르(modified Picard) 구조를 쓴다. 정확도가 아니라 보존성을 이산화 단계에서 확보한 사례로 자주 인용된다.

6. 피카르인가 뉴턴인가[편집]

수정 피카르는 $K$ 와 $C$ 를 직전 반복값으로 얼려 쓴다. 즉 계수의 $h$ 의존성을 야코비안에 넣지 않는다.

• 수정 피카르 — 행렬이 항상 대칭 양정치라 켤레기울기법을 그대로 태울 수 있다. 수렴은 1차(선형). 구현이 짧고 잘 안 죽는다.
• 뉴턴-랩슨법 — $\partial K/\partial h$, $\partial C/\partial h$ 항까지 자코비안 행렬에 넣는다. 2차 수렴이지만 행렬이 비대칭이 되고, 초기 추정이 나쁘면 발산한다.

파니코니와 푸티(1994)의 비교가 자주 인용되는 결론을 준다: 완만한 문제에서는 피카르가 반복당 비용이 싸서 이기고, 극도로 건조한 초기조건이나 급한 전선에서는 피카르가 아예 수렴을 못 해 뉴턴만 답을 낸다. 그래서 실무 코드는 피카르로 시작해 반복이 밀리면 뉴턴으로 갈아타는 하이브리드를 쓴다. 뉴턴을 쓸 때는 선탐색이나 신뢰영역이 사실상 필수인데, 이유가 재미있다 — 반 헤누흐텐-무알렘 조합은 $n<2$ 인 세립 토양에서 $S_e\to1$ 일 때 $\mathrm{d}K/\mathrm{d}S_e$ 가 발산한다. 점토처럼 $n\approx1.1$ 인 토양이 정확히 그 경우라, 포화 근처에서 뉴턴 스텝이 그대로 폭주한다. 포화 직전에 작은 공기침입값을 강제로 끼워 넣어 곡선을 자르는 보정(Vogel 등, 2001)이 널리 쓰이는 것은 물리보다 수렴 때문이다.

시간 간격은 거의 항상 반복 횟수로 제어한다. “3회 이내 수렴이면 $\Delta t$ 를 1.3배, 10회 초과면 절반으로 줄이고 스텝 재시도”류의 휴리스틱이며, 상한은 격자 하나를 전선이 통과하는 시간으로 잡는다. 그리고 코드가 매 스텝 반드시 출력해야 하는 값이 질량균형비(누적 저장 변화 ÷ 누적 순유입)다. 이게 1에서 벗어나면 계산이 예쁘게 그려져도 답이 아니다.

수렴 자체를 보장하려는 방향도 있다. 야코비안 대신 충분히 큰 상수 $L$ 을 쓰는 L-도식은 초기값과 $\Delta t$ 에 무관하게 수축사상임이 증명되어 있지만, 수렴이 1차이고 느리다. 견고함과 속도의 교환이라는 구도가 여기서도 반복된다.

7. 경계조건이 스스로 종류를 바꾼다[편집]

리처즈 방정식 구현에서 초보자가 가장 많이 걸리는 곳은 솔버가 아니라 경계조건이다. 지표면 조건이 해에 따라 종류가 바뀌기 때문이다.

비가 강도 $i$ 로 내린다고 하자. 처음에는 내리는 대로 다 스미므로 지표는 규정 플럭스(노이만) 조건이다.

$$-K(h)\left(\frac{\partial h}{\partial z}+1\right)\Big|_{z=0}=i$$

그런데 토양이 그 유량을 못 받아들이게 되는 순간이 온다. 그러면 물이 고이고 지표 압력수두가 고인 깊이로 정해지므로 조건이 규정 수두(디리클레) $h=h_0\approx0$ 으로 전환된다. 반대로 증발 상황에서는 대기와 평형인 수두 하한($h\approx-10^4\sim-10^5\,\mathrm{cm}$)에 걸리는 순간 다시 전환된다. 즉

$$q\le q_{\max},\qquad h\le h_0,\qquad (q-q_{\max})(h-h_0)=0$$

라는 상보 조건이며, 수학적으로는 변분부등식이다. 코드는 매 반복마다 각 지표 절점의 조건 종류를 판정하는데, 이 판정이 반복 사이에서 왔다 갔다 하면(채터링) 수렴 판정이 영원히 안 붙는다. 조건이 바뀐 절점 수를 세어 강제로 고정하거나, 스텝을 줄여 다시 시도하는 처리가 필요하다.

같은 구조가 하부·측면에도 나온다.

• 자유배수(free drainage) — 바닥에서 $\partial h/\partial z=0$ 을 걸면 $q=-K(h)$, 즉 중력만으로 빠져나간다. 지하수면이 깊을 때의 표준 하부 조건.
• 지하수면 — $h=0$ 인 디리클레.
• 침윤면(seepage face) — 제방·사면의 하류측 비탈. 물이 나오는 면에서는 $h=0$, 안 나오는 면에서는 $q=0$ 인데 그 경계가 어디인지가 미지수다. 지표 조건과 똑같은 상보 구조이며, 제방 침투 해석이 유난히 잘 안 붙는 이유다.

8. 히스테리시스[편집]

같은 압력수두인데 적시는 중이냐 마르는 중이냐에 따라 함수량이 다르다. 마르는 곡선(배수)이 적시는 곡선(흡습)보다 위에 있고, 그 사이를 오가는 부분 곡선을 주사곡선(scanning curve)이라 한다. 원인은 셋이다.

• 잉크병 효과 — 목이 좁고 배가 부른 공극은, 배수될 때는 좁은 목의 모관력이 버티고 흡습될 때는 넓은 배가 먼저 채워지지 않는다. 순수한 기하학적 비대칭.
• 접촉각 히스테리시스 — 전진 접촉각이 후퇴 접촉각보다 크다.
• 갇힌 공기 — 흡습 시 공극에 공기가 갇혀 흡습 곡선의 최대 함수량이 배수 곡선의 $\theta_s$ 에 못 미친다. 보통 0.85~0.95배.

실무적으로 중요한 사실 하나: $K(\theta)$ 의 히스테리시스는 $\theta(h)$ 에 비해 훨씬 약하다. 즉 전도도는 함수량만 알면 거의 정해지고, 압력수두와 함수량을 잇는 관계만 이력에 의존한다. 그래서 히스테리시스를 넣을 때는 보유곡선만 이력 모형으로 바꾸고 $K(\theta)$ 는 단일 곡선을 쓰는 것이 표준이다.

구현은 두 주곡선 사이에서 주사곡선을 스케일링으로 만드는 방식(Kool–Parker)이 가장 흔하다. 문제는 각 셀이 자신의 반전 이력을 상태변수로 들고 있어야 한다는 것 — 방정식이 더 이상 $(h,\theta)$ 만의 함수가 아니게 되어 메모리와 분기가 늘고, 반복 중에 반전점이 바뀌면 수렴이 흔들린다. 그래서 히스테리시스가 중요한 것을 다들 알면서도 대부분의 모의가 배수 곡선 하나로 끝난다. 강우-건조가 반복되는 연속 모의에서 이 생략은 첨두 함수량을 조용히 틀리게 만든다.

9. 리처즈가 못 하는 것[편집]

방정식 자체의 가정에서 오는 한계도 있다.

• 공기상을 무시했다. 기체 압력을 대기압 상수로 뒀는데, 배기가 막힌 층상 토양에서는 갇힌 공기가 압축되어 침투를 억제한다. 물-공기 2상 방정식을 제대로 풀어야 잡힌다.
• 중력 핑거링을 못 낸다. 마른 모래에 물을 부으면 전선이 평평하게 내려가지 않고 손가락처럼 갈라지는데, 리처즈 방정식은 균질 매질에서 이 불안정을 만들지 못한다. 포물형이라 해가 매끄럽게 뭉개지기 때문. 하사니자데-그레이의 동적 모관압 항($p_c$ 가 $\partial\theta/\partial t$ 에도 의존)을 추가해야 재현된다.
• 거대공극·균열. 뿌리 구멍과 건조 균열을 타고 물이 전선을 추월한다. 매질을 둘로 쪼개는 이중공극 모형으로 우회한다.
• REV 가정. 애초에 다르시 법칙이 성립할 대표체적요소가 있어야 한다.

10. 관련 문서[편집]

• 다르시 법칙 · 다공성 매질 유동 · 그린-암프트 모형
• 토양수분특성곡선 · 수리전도도 · HYDRUS
• 뉴턴-랩슨법 · 자코비안 행렬 · 준-뉴턴법
• 대류-확산 방정식 · 강성 방정식 · 유한체적법
• SWMM · 검증 및 확인 · 불확실성 정량화
• 표면장력 · 다상유동

11. Footnotes[편집]

1. 이름 문제가 하나 있다. 이 방정식은 1922년 루이스 프라이 리처드슨(Lewis Fry Richardson)의 『Weather Prediction by Numerical Process』에 이미 등장하고, 1931년 로렌조 리처즈(Lorenzo A. Richards)가 독립적으로 유도해 실험과 함께 발표하면서 그 이름이 붙었다. 최근 문헌은 리처드슨-리처즈 방정식으로 병기하는 쪽이 늘고 있다. 참고로 리처드슨은 전산유체역학에서도 “인간 6만 4천 명 일기예보 공장”으로 등장하는 그 사람이다. 이 양반 한 명이 수치기상예보와 토양물리 방정식을 동시에 선점했다. ↩

2. 그래서 반 헤누흐텐 매개변수를 현장에서 재는 대신 토성(모래·실트·점토 비율)과 용적밀도에서 회귀로 뽑아 쓰는 토양전달함수(pedotransfer function)가 널리 쓰인다. 흙 시료를 체로 쳐서 얻은 세 숫자로 $\alpha$ 와 $n$ 을 추정하고, 그걸로 $K$ 를 몇 자릿수 범위에서 예측하겠다는 뜻이다. 잘 될 리가 없고 실제로 잘 안 되지만, 대안이 압력판 실험 몇 주라서 다들 쓴다. ↩

3. 그래서 리처즈 방정식 논문의 벤치마크 목록은 20년째 거의 같다. Celia의 건조 모래 침투, 폴만스 배수 문제, 트레이시의 해석해 정도. 새 도식을 만들면 이 셋으로 두들겨 맞고, 대부분 “우리 도식은 Celia 문제에서 CPU 시간이 40% 짧습니다”라는 문장으로 끝난다. 그 40%가 다른 토양에서 유지된다는 보장은 아무도 안 한다. ↩