\underline{寒冷氷床の離散モデル化目次 1}

寒冷氷床の離散モデル化

村田泰洋

2002 年 08 月 20 日

目次

1 座標系 2 1.1 鉛直座標変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 格子点の設定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 支配方程式 5 2.1 速度場 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 氷厚の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 熱力学の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

\underline{寒冷氷床の離散モデル化  1 座標系  2}

1 座標系

ここでは計算する際に必要な鉛直方向の座標変換, および空間格子点の切り方を記す.

1.1 鉛直座標変換

数値的に計算するのに鉛直座標をその場所の氷厚 H(x, t) = h(x, t) − b(x, t))) で引き延ばされた無次元座標 ζ

ζ = \frac{z − b(x, t)}{H(x, t)}

に変換する (Jenssen, 1977). ここで h(x, t) は表面高度, b(x, t) は基盤高度である)

この座標系では氷床上部境界が常に格子点上にあるので上部境界の数値誤差を減らすことができ, また氷厚に従って座標が貼り変えられているので氷床の縁ではより正確な計算を行う事ができる. ^{1}.

ζ 座標を用いることで微分, 積分は以下のように変換される. 関数 f は変数 x, z(ζ), t の任意の関数とし, r は x, t のいずれかである.

(\frac{∂f}{∂r} )_{z} = (\frac{∂f}{∂r} )_{ζ} − \frac{1}{H} (\frac{∂f}{∂ζ} ) (\frac{∂b}{∂r} + ζ \frac{∂H}{∂r} )
\frac{∂f}{∂z} = \frac{1}{H} \frac{∂f}{∂ζ}
\frac{∂^{2}f}{∂z^{2}} = \frac{1}{H^{2}} \frac{∂^{2}f}{∂ζ^{2}}
∫ f(z) dz = H ∫ f(ζ) dζ

これによって支配方程式は以下のようになる.

\frac{∂H}{∂t} = −\frac{∂F}{∂x} + \frac{M_{s}}{ρ}  (1.1)

\frac{∂T}{∂t} − \frac{1}{H} \frac{∂T}{∂ζ} (\frac{∂b}{∂t} + ζ \frac{∂H}{∂t} )

= − (u (\frac{∂T}{∂x} − \frac{1}{H} \frac{∂T}{∂ζ} (\frac{∂b}{∂x} + ζ \frac{∂H}{∂x} )) + \frac{w}{H} \frac{∂T}{∂ζ} )

1しかし実際に計算した際に熱力学の式を解く計算において氷床の縁でよくわからない数値的なゴミが発生した.

\overline{main.tex 2002 年 08 月 20 日 (村田泰洋)}

\underline{寒冷氷床の離散モデル化  1 座標系  3}
+\frac{k}{ρc_{p}H^{2}} \frac{∂^{2}T}{∂ζ^{2}} + \frac{2}{ρc_{p}} A(T)σ_{zx}^{n+1}  (1.2)
F = −2 (\frac{ρg}{η} )^{n} |\frac{∂h}{∂x}|^{n−1} \frac{∂h}{∂x}H^{n+2} ∫_{0}^{1} dζ ∫_{0}^{ζ} A(T)(1 − ζ′)^{n} dζ′  (1.3)
u = −2 (\frac{ρg}{η} )^{n} |\frac{∂h}{∂x}|^{n−1} \frac{∂h}{∂x}H^{n+1} ∫_{0}^{ζ} A(T)(1 − ζ′)^{n} dζ′  (1.4)
w = −H ∫_{0}^{ζ} \frac{∂u}{∂x} dζ + \frac{∂H}{∂x} (uζ − ∫_{0}^{ζ} u dζ) + \frac{∂b}{∂x} (u(ζ) − u(b))))  (1.5)
σ_{xz} = −ρg \frac{∂h}{∂x}h(1 − ζ)  (1.6)
A(T) = A_{0} exp (−\frac{Q}{RT} )  (1.7)

1.2 格子点の設定

格子点は図 (1.1)) のように温度, 氷厚を整数格子点に, 速度, フラックスを半整数格子点に取る. 格子点以外での値は線形補間によって表す^{2}. 水平格子水平方向の格子数は N_{i} とし, 格子間隔を等間隔 Δx とする. 鉛直格子鉛直方向の格子点は整数格子点と半整数格子点を定義し, それぞれの格子数 N_{j} 個とし, 格子間隔は非一様とする. 格子間隔は整数格子点, 半整数格子点を用いて

Δζ_{j} = ζ_{j+ \frac{1}{2}} − ζ_{j− \frac{1}{2}}  (j = 1, 2, ⋯ N_{j})  (1.8)
Δζ_{j+ \frac{1}{2}} = ζ_{j+1} − ζ_{j}  (j = 1, 2, ⋯ N_{j−1})  (1.9)

と定義する. ここで ζ_{\frac{1}{2}} , ζ_{N_{j}+ \frac{1}{2}} はそれぞれ下部境界, 上部境界とする.

2例えば h_{i+1/2} = \frac{h_{i+1}+h_{i}}{2} \overline{main.tex 2002 年 08 月 20 日 (村田泰洋)}

図 1.1: 変数の格子点位置. 上の図が空間 1 次元の変数で, 下の図が空間 2 次元.

\overline{main.tex 2002 年 08 月 20 日 (村田泰洋)}

\underline{寒冷氷床の離散モデル化  2 支配方程式  5}

2 支配方程式

2.1 速度場

速度場の式 (1.4), (1.5) を水平方向には前方差分によって離散化すると

u_{i+ \frac{1}{2}} = −2 (\frac{ρg}{τ} )^{n} |\frac{h_{i+1} − h_{i}}{Δx} |^{n−1} \frac{h_{i+1} − h_{i}}{Δx} (H_{i+ \frac{1}{2}} )^{n+1}
× ∫_{0}^{ζ} A(T)(1 − ζ′)^{n} dζ′  (2.10)
w_{i} = −H_{i} ∫_{0}^{ζ} \frac{u_{i+ \frac{1}{2}} − u_{i− \frac{1}{2}}}{Δx} dζ′

+\frac{h_{i+ \frac{1}{2}} − h_{i− \frac{1}{2}}}{Δx} (u_{i}ζ + ∫_{0}^{ζ} u_{i} dζ′)

+\frac{b_{i+ \frac{1}{2}} − b_{i− \frac{1}{2}}}{Δx} (u_{i} − u_{0}))  (2.11)

となる. 鉛直方向の離散化は台形公式を用いて以下のように表す^{3} .

∫_{0}^{ζ} I dζ′ = \frac{ζ′_{1}}{2} I_{1} + ∑_{j′=2}^{j} \frac{Δζ_{j′+ \frac{1}{2}}}{2} (I_{j′} + I_{j′−1}) (j = 2, 3, ⋯ N_{j}) (2.12)

ここで I は u もしくは A(1 − ζ)^{n} である.

2.2 氷厚の式

氷厚の式は空間方向に前方差分, 時間方向に陽解法を用いて離散化する.

\frac{H_{i}^{k+1} − H_{i}^{k}}{Δt} = −\frac{F_{i+ \frac{1}{2}}^{k} − F_{i− \frac{1}{2}}^{k}}{Δx} + \frac{M_{i}^{k}}{ρ}  (2.13)

ここで F_{i+ \frac{1}{2}} は速度場の式の鉛直離散化を施した時と同様に台形公式を用いて求める.

3ζ = 0, ζ = 1 のときの値が記されていないので追記すべし.

\overline{main.tex 2002 年 08 月 20 日 (村田泰洋)}

\underline{寒冷氷床の離散モデル化  2 支配方程式  6}

F_{i+ \frac{1}{2}} = \frac{ζ_{1}}{2} I_{1} + ∑_{j′=2}^{N_{j}} \frac{Δζ_{j′+ \frac{1}{2}}}{2} (I_{j′} + I_{j′−1}) + \frac{ζ_{N_{j}+ \frac{1}{2}} − ζ_{N_{j}}}{2} (I_{N_{j}+ \frac{1}{2}} + I_{N_{j}}) (2.14)

2.3 熱力学の式

熱力学の式は Huybrechts 1990, saito 2002 で用いられたスキームに従う^{4}.

空間方向に関して鉛直方向に中心差分, 水平方向には風上差分を用いて離散化する.

\frac{∂T}{∂t} − \frac{V_{i,j}}{H_{i}} (\frac{T_{i,j+1} − T_{i,j−1}}{2Δζ_{j}} )

− \frac{k}{ρc_{p}} \frac{1}{(H_{i})^{2}} \frac{1}{Δζ_{j}} (\frac{T_{i,j+1} − T_{i,j}}{Δζ_{j+ \frac{1}{2}}} − \frac{T_{i,j} − T_{i,j−1}}{Δζ_{j− \frac{1}{2}}} )

= −u_{i,j} \frac{T_{i,j} − T_{i−1,j}}{Δx} + \frac{1}{ρc_{p}} A(T_{i,j})(σ_{xz})_{i,j}^{n+1}  (2.15)

左辺第 2 項の V_{i,j} は

V_{i,j} = \frac{∂b}{∂t} + ζ_{j} \frac{∂h}{∂t} + u_{i,j} (ζ_{j} \frac{H_{i+1} − H_{i−1}}{2Δx} + \frac{b_{i+1} − b_{i−1}}{2Δx} ) − w_{i,j}  (2.16)

とする. また右辺第 2 項の (σ_{xz})_{i,j} は

(σ_{xz})_{i,j} = −ρgh_{i} \frac{h_{i+1} − h_{i−1}}{2δx} (1 − ζ_{j})  (2.17)

となる.

時間方向には鉛直微分に関しての陰解法を用いて離散化する.

\frac{T^{k+1} − T^{k}}{Δt} − \frac{V^{k}}{h^{k}} (\frac{∂T}{∂ζ} )^{k+1} − \frac{k}{ρc_{p}} \frac{1}{h^{k}} (\frac{∂^{2}T}{∂ζ^{2}} )^{k+1}

4このスキームは

• 鉛直微分は中心差分, 水平微分は風上差分できていない.

• 時間微分は温度の鉛直微分に関して陰解法を用いる. このスキームはフォン・ノイマンの安定条件を満たしてはいる (証明は省略) が, 物理的制約条件を満たしているのかがどこにも示されていないため, 真に有用なスキームであるかは確認 \overline{main.tex 2002 年 08 月 20 日 (村田泰洋)}

= −u^{k} (\frac{∂T}{∂x} )^{k} + \frac{1}{ρc_{p}} A(T^{k})σ_{xz}^{k}^{n+1} (2.18)
V^{k} = ζ \frac{h^{k+1} − h^{k}}{δt} + uζ (\frac{∂h}{∂x})^{k} − w^{k}  (2.19)

以上を整理すると

T_{i,j+1}^{k+1}\{\frac{−V_{i,j}^{k+1}}{h_{i}^{k}} \frac{Δt}{2Δζ_{j}} − \frac{k}{ρc_{p}} \frac{Δt}{(h_{i}^{k})^{2}} \frac{1}{Δζ_{j}} \frac{1}{Δ_{j+ \frac{1}{2}}} \}

+ T_{i,j}^{k+1}\{1 + Δt\frac{k}{ρc_{p}} \frac{1}{(h_{i}^{k})^{2}} \frac{1}{Δζ_{j}} (\frac{1}{Δζ_{j+ \frac{1}{2}}} + \frac{1}{Δζ_{j− \frac{1}{2}}} )\}

+ T_{i,j−1}^{k+1}\{\frac{V_{i,j}^{k}}{h_{i}^{k}} \frac{Δt}{2Δζ_{j}} − \frac{k}{ρc_{p}} \frac{Δt}{(h_{i}^{k})^{2}} \frac{1}{Δζ_{j}} \frac{1}{Δζ_{j− \frac{1}{2}}} \}

= T_{i,j}^{k} (1 − \frac{Δt}{Δx}u_{i,j}^{k}) + T_{i−1,j}^{k} (\frac{Δt}{Δx}u_{i,j}^{k}) + \frac{1}{ρc_{p}} A(T)σ_{xz}^{n+1} (2.20)
V_{i,j}^{k} = ζ_{j} \frac{h_{i}^{k+1} − h_{i}^{k}}{Δt} + u_{i,j}^{k}ζ_{j} \frac{h_{i+1}^{k} − h_{i−1}^{k}}{2Δx} − w_{i,j}^{k}  (2.21)
(σ_{xz})_{i,j}^{k} = −ρgh_{i}^{k} \frac{h_{i+1}^{k} − h_{i−1}^{k}}{2δx} (1 − ζ_{j})  (2.22)

となる.

参考文献

• Jenssen, D., 1977. A three-dimensional polar ice-sheet model. J. Glaciol, 18, 373 - 389.

• Huybrechts,P.1990. A 3-D model for the Anterctic ice sheet: a sensitivity study on the glacial-interglacial contrast. Climate Dynamics 5, 79–92

• F.saito. 2002. Development of a Three Dimensional Ice Sheet Model for Numerical Studies of Antarctic and Greenland Ice Sheet CCSR Report 15, pp135