원시방정식

편집 역사 토론
유체역학 전산유체역학 계산물리 마지막 수정: 2026-08-04 04:51:20

1. 개요[편집]

원시방정식(primitive equations)은 회전하는 구면 위의 얇은 대기층에 대해 나비에-스토크스 방정식을 정수압 근사와 얕은 대기 근사로 줄여 얻은 연립 편미분방정식계로, 수평 운동량 방정식·질량 보존·열역학 제1법칙·이상기체 상태방정식·수증기 보존을 한 벌로 묶은 것이다. 지구 대기를 도메인으로 쓰는 거의 모든 모델의 역학코어(dynamical core)가 이 방정식, 혹은 정수압 가정만 풀어 준 비정수압 확장판을 시간 적분한다.

이름의 “원시”는 원시적이라는 뜻이 아니라 원변수(primitive variable) — 속도·온도·기압·비습 그 자체 — 로 쓰였다는 뜻이다.1 와도와 유선함수 같은 유도 변수로 쓰인 준지균 모형이나 얕은 물 방정식 계열의 필터 모형과 대비되는 이름이다.

이 문서는 지배방정식과 그 이산화만 다룬다. 이걸 업무 예보 시스템으로 조립하는 이야기는 수치기상예보, 초기장을 만드는 이야기는 자료동화, 확률 예보로 확장하는 이야기는 앙상블 예보에 위임한다.

2. 무엇을 버려서 얻은 방정식인가[편집]

출발점은 회전계에서 쓴 완전한 압축성 유체 방정식이다. 여기에 세 개의 칼이 들어간다.

얕은 대기 근사(shallow atmosphere). 대기의 두께는 10 km 남짓, 지구 반지름은 6371 km다. 비가 100배가 넘으니 반지름 rr을 상수 aa로 고정하고, 그로부터 파생되는 계량항(metric term)들을 함께 버린다.

전통적 근사(traditional approximation). 코리올리 항 2Ω×v2\boldsymbol{\Omega}\times\mathbf{v} 중 연직속도 ww와 곱해지는 2Ωcos⁡φ2\Omega\cos\varphi 성분을 버리고 f=2Ωsin⁡φf = 2\Omega\sin\varphi만 남긴다(Ω=7.292×10−5 s−1\Omega = 7.292\times10^{-5}\ \mathrm{s^{-1}}). 이건 얕은 대기 근사와 세트로 버려야 각운동량·에너지 보존이 깨지지 않는다는 것이 필립스(1966)가 정리한 요점이다.2

정수압 근사(hydrostatic approximation). 연직 운동방정식을 통째로 진단 관계식으로 대체한다.

∂p∂z=−ρg\frac{\partial p}{\partial z} = -\rho g

연직 가속도 Dw/DtDw/Dt가 중력 대비 (H/L)2\left(H/L\right)^2 크기라는 척도 해석이 근거다. 종관 규모(L∼1000L\sim 1000 km, H∼10H\sim 10 km)에서 이 값은 10−410^{-4} 수준이라 마음 편히 버릴 수 있다. 대가는 뒤에서 청구된다.

3. 방정식 한 벌[편집]

기압을 연직좌표로 쓰면(ω=Dp/Dt\omega = Dp/Dt, Φ=gz\Phi = gz는 지오퍼텐셜) 원시방정식은 다음의 다섯 줄로 정리된다.

DvhDt+f k^×vh=−∇pΦ+F\frac{D\mathbf{v}_h}{Dt} + f\,\hat{\mathbf{k}}\times\mathbf{v}_h = -\nabla_p \Phi + \mathbf{F} ∂Φ∂p=−RTp,∇p⋅vh+∂ω∂p=0\frac{\partial \Phi}{\partial p} = -\frac{RT}{p}, \qquad \nabla_p\cdot\mathbf{v}_h + \frac{\partial \omega}{\partial p} = 0 DTDt−κTp ω=Qcp,DqDt=Sq,κ=Rcp≈0.286\frac{DT}{Dt} - \frac{\kappa T}{p}\,\omega = \frac{Q}{c_p}, \qquad \frac{Dq}{Dt} = S_q, \qquad \kappa = \frac{R}{c_p} \approx 0.286

여기서 D/Dt=∂t+vh⋅∇p+ω ∂pD/Dt = \partial_t + \mathbf{v}_h\cdot\nabla_p + \omega\,\partial_p다. 눈여겨볼 점은 연직속도가 예단 변수가 아니라는 것이다. ω\omega는 연속방정식을 연직으로 적분해서 진단하고, 지오퍼텐셜은 정수압 관계를 적분해서 진단한다. 연직 방향에는 시간 미분이 남아 있지 않다.

로스비수 Ro=U/fLRo = U/fL가 작은 종관 규모에서는 좌변의 코리올리 항과 우변의 기압경도력이 거의 상쇄되어 지균 균형에 가까운 상태가 된다. 실제 대기는 이 균형에서 몇 퍼센트만 벗어나 있고, 날씨는 바로 그 몇 퍼센트에서 나온다. 이 사실이 이산화 전략 전체를 좌우한다 — 큰 두 항의 차이를 보는 계산이므로, 격자가 두 항을 서로 다른 정확도로 계산하면 결과가 망가진다. 코리올리 항과 기압경도항을 같은 격자점 배치에서 일관되게 다루기 위해 엇갈림 격자(아라카와 C-격자 등)를 쓰는 이유가 여기 있다. 회전계 자체의 일반론은 회전 좌표계 기법을 참고.

4. 연직좌표계 — 지형을 어디에 숨길 것인가[편집]

기압좌표는 깔끔하지만 치명적 결함이 있다. 산이 있으면 하부 등압면이 땅속으로 들어간다. 그래서 실제 모델은 지형을 좌표계 안으로 흡수한다.

좌표정의성격
시그마 σ\sigmaσ=p/ps\sigma = p/p_s지형 추종, 하부 경계가 항상 σ=1\sigma=1
하이브리드 η\etapk=Ak+Bk psp_k = A_k + B_k\,p_s하부는 시그마, 상부는 기압
고도 zz / 지형추종 z∗z^*기하 고도 기반비정수압 모델에서 표준
등온위 θ\theta온위 기반단열 유동이 좌표면 위를 미끄러짐

필립스(1957)의 시그마 좌표는 하부 경계조건을 σ˙=0\dot\sigma = 0으로 만들어 버리는 우아함이 있지만, 기압경도력이 두 항으로 쪼개진다.

−∇pΦ  =  −∇σΦ−RT ∇σln⁡ps-\nabla_p \Phi \;=\; -\nabla_\sigma \Phi - RT\,\nabla_\sigma \ln p_s

가파른 지형 위에서는 이 두 항이 각각 거대하고 거의 상쇄된다. 그러면 큰 수에서 큰 수를 빼는 전형적인 상황이 되어, 각 항의 절단오차가 잔차보다 커지는 사고가 난다. 티베트 고원 위에서 가짜 바람이 부는 유명한 증상이 이것이다. 하이브리드 좌표(AkA_k, BkB_k를 잡아 하부만 지형 추종, 상부는 순수 기압면)는 이 오차를 성층권까지 끌고 가지 않으려는 처방이고, 시몬스-버리지(1981) 이후 전지구 모델의 사실상 표준이 되었다. 요즘 업무 모델은 50~140층 정도를 쓴다.

5. 음파를 지운 대가[편집]

정수압 근사의 진짜 목적은 지형 처리도 물리 단순화도 아니고 연직 전파 음파의 제거다. 음속은 340 m/s인데 기상학적으로 의미 있는 신호는 수십 m/s다. 음파를 그대로 두면 CFL 조건이 아무 정보도 나르지 않는 파에 의해 결정된다.

다만 정수압계가 음파를 전부 지우지는 않는다. 수평으로 달리는 램 파(Lamb wave)는 음속에 가까운 속도로 살아남고, 중력파도 그대로 남는다. 음파를 완전히 죽이고 싶으면 비탄성 근사(anelastic)나 유사비압축 근사로 가야 하는데, 그 순간 매 스텝 타원형 포아송 방정식을 풀어야 해서 비용 구조가 바뀐다. 세상에 공짜는 없다.

더 비싼 대가는 불균형에 대한 취약성이다. 원시방정식은 준지균 모형과 달리 중력파를 그대로 품고 있어서, 초기장이 지균 균형에서 조금만 벗어나 있어도 대진폭 중력파가 터진다. 리처드슨의 1922년 손계산이 6시간 만에 지상기압 145 hPa 변화를 뱉은 원인이 정확히 이것이다. 비선형 정규모드 초기화(Machenhauer, 1977)와 그 후신인 현대 자료동화의 균형 제약은 이 문제에 대한 답이다.

6. 정수압이 깨지는 지점 — 수 km 문턱[편집]

정수압 근사의 유효 조건은 운동의 수평 규모가 연직 규모보다 훨씬 크다는 것이다. 깊은 대류 셀은 폭 1~10 km, 깊이 10 km라서 이 조건이 성립하지 않는다. 실무 감각으로 정리하면 대략 이렇다.

  • 격자 Δx≳10\Delta x \gtrsim 10 km: 정수압으로 충분. 대류는 파라미터화한다.
  • Δx≈3–10\Delta x \approx 3\text{–}10 km: 회색지대(grey zone). 대류가 부분적으로 해상되는데 파라미터화도 아직 필요한, 아무도 좋아하지 않는 구간.
  • Δx≲3\Delta x \lesssim 3 km: 비정수압 필수. 대류 허용(convection-permitting) 영역.

전지구 모델이 10 km 아래로 내려온 지금은 애초에 비정수압 코어로 짓는 것이 국룰이 되었다. 그래서 “요즘도 원시방정식을 쓰냐”는 질문의 답은 쓴다 — 연직 운동방정식을 완전판으로 되돌리고 나머지 구조(구면 좌표, 하이브리드 연직좌표, 열역학·수증기 방정식, 균형 문제)는 그대로 물려받기 때문이다.

7. 실제로 푸는 방법[편집]

스펙트럴 변환법. 수평 방향을 구면조화 함수로 전개하면 미분이 대수 연산이 되고 극점 수렴 문제가 사라진다. 다만 비선형 항은 격자에서 계산해야 해서 매 스텝 스펙트럴↔격자 변환이 필요하다. 경도 방향은 고속 푸리에 변환으로 싸게 처리되지만 위도 방향 르장드르 변환은 절단파수 NN에 대해 O(N3)O(N^3)이라, 해상도를 올릴수록 이쪽이 비용을 지배한다. 고속 르장드르 변환이 나오기 전까지 이것이 스펙트럴 모델의 해상도 천장이었다. 자세한 배경은 스펙트럴 방법 참고. 대안인 격자점 모델은 유한체적법이나 분광요소법을 정20면체·육면체구 격자에 얹어 통신 부하를 줄인다.

세미-라그랑주 기법. 이류를 도착점에서 역궤적을 추적해 출발점 값을 보간하는 방식으로 풀면 이류 CFL 제약이 사라진다. 제트기류가 100 m/s로 불어도 시간 간격을 유지할 수 있다는 뜻이다. 단점은 질량 보존이 자동이 아니라는 것과, 보간 차수가 곧 수치 소산과 분산 특성을 결정한다는 것. 이걸 세미-임플리시트와 결합한 SISL 조합이 스펙트럴 전지구 모델의 오랜 국룰이다.

시간 분할(split-explicit). 클렘프-빌헬름슨(1978)이 제안한 반대 방향의 해법이다. 느린 모드(이류·물리)는 큰 시간 간격 Δt\Delta t로, 음파와 중력파에 관련된 항만 작은 부분 스텝 Δτ\Delta\tau로 여러 번 반복해 적분한다. 큰 행렬을 푸는 대신 싼 스텝을 4~6번 더 밟는 셈이라, 타원형 solver를 피하고 싶은 비정수압 지역 모델(WRF, COSMO, MPAS 계열)이 즐겨 쓴다.3

셋 중 무엇을 고르든 마지막에 남는 문제는 같다. 보존이다. 스펙트럴 변환은 질량·에너지를 자동으로 보존하지 않고, 세미-라그랑주는 더 심하며, 시간 분할은 부분 스텝 사이에서 균형이 어긋난다. 그래서 실제 역학코어에는 항상 질량 고정자(mass fixer)와 에너지 보정 항이 붙어 있다. 예보 3일 차의 오차보다 100년 기후적분에서의 질량 누출이 더 무서운 이유다.4

8. 관련 문서[편집]

9. Footnotes[편집]

  1. 그래서 학부생이 “원시방정식”이라는 번역어를 처음 보고 “그럼 문명 방정식도 있나요”라고 묻는 것은 통과의례에 가깝다. 영어로도 primitive라는 단어 선택이 좋진 않았다는 데 대체로 동의하는 분위기다. ↩

  2. 이걸 어긴 근사를 쓰면 방정식계가 각운동량을 보존하지 않는다. 하루 이틀 예보에서는 티가 안 나지만 수십 년 기후적분을 돌리면 대기가 서서히 회전 상태를 잃어버린다. 근사는 세트로 사야 한다. ↩

  3. 음파 부분 스텝은 보통 전진-후진(forward-backward) 방식으로 짠다. 완전 양해적으로 두면 부분 스텝의 CFL이 다시 음속에 묶여서, 애초에 분할한 의미가 없어지기 때문이다. ↩

  4. 100년 적분에서 수증기 총량이 슬금슬금 새는 모델은 어느 시점부터 지구를 사막으로 만든다. 그래서 기후 모델러들은 “예쁜 예보”보다 “새지 않는 코어”를 훨씬 중요하게 본다. ↩