\underline{寒冷氷床モデルの定式化  1 はじめに  1}

寒冷氷床モデルの定式化

村田泰洋

2001 年 07 月 11 日

目次

1 はじめに 1

2 寒冷氷床モデルの定式化 2 2.1 考える系 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 氷床力学の基礎方程式系 . . . . . . . . . . . . . . . . . . . . . . . . 2 2.3 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.3.1 運動学的境界条件 . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3.2 力学的境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.3.3 熱学的境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.4 直線直行座標での表示 . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.5 浅氷近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.5.1 運動量保存則 . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.5.2 エネルギー保存則 . . . . . . . . . . . . . . . . . . . . . . . . 7 2.5.3 構成方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.5.4 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1 はじめに

私は氷床外部の気候条件によって氷床がどのように振舞うか, ということに興味を持っている.

\underline{寒冷氷床モデルの定式化  2 寒冷氷床モデルの定式化  2}

具体的な話としては例えば, 平均気温の変化が氷床の質量収支にどのように影響を与えているかなどがそれに当たる. このような具体的な問題に対して地球上に存在する各々の氷河・氷床での応答実験は行われている. しかし現実に存在する氷床の性質を理解するためには簡単な枠組で考えることが有効である (かもしれない). そこで仮想的な氷床と気候条件を設定してそこでの氷床の基本的な振舞い (例えば, 氷厚・温度・流速など) について理解する.

今回の M ゼミでは氷床流動モデルの中で最も基本的なモデルである寒冷氷床について注目し, この氷床の定式化を行った.

今後の計画としては

1. 定式化によって得られる方程式系を離散化して数値計算コードを作る

2. 用いたコードの結果と EISMINT が行ったベンチマーク実験と比較する

3. さまざまな構成方程式 (例えば異方性のある氷の構成方程式や CO_{2} の構成方程式) を用いたの氷床の振舞いについて計算する.

の順に行っていく.

2 寒冷氷床モデルの定式化

2.1 考える系

図 2.1 のような水平方向に無限に広がる空間に年間涵養量 M_{s} で成長する氷体 (ice body) を考える. この氷体は非ニュートン・非圧縮・熱伝導流体 (Hutter 1983) として扱う. 氷体の温度は十分に低いものとして融解は起こらないとする. 氷体の水平厚さ L は氷厚 H に比べて十分大きいものとする. また氷体の重みによる基盤地形の変形はないと考える. 水平方向の座標を x, y, 鉛直方向を z とする.

2.2 氷床力学の基礎方程式系

基礎方程式系は次の質量, 運動量, エネルギーの保存則からなる.

divv = 0  (2.1)
0 = divσ + ρg  (2.2)
ρ\frac{dc_{p}T}{dt} = −div(kgradT) + tr(σ · D)  (2.3)

図 2.1: 系の設定

氷体の流動が非常にゆっくりなので運動量保存則右辺は 0 とする. 記号の意味はそれぞれ次の通りである.

\begin{matrix}v & 速度ベクトル = (u, v, w) \\ ρ & 氷の密度 \\ σ & 応力テンソル \\ g & 重力加速度 \\ T & 氷温 \\ D & 歪み速度テンソル \\ k & 熱伝導度 \\ c_{p} & 氷の比熱 \end{matrix}

構成方程式は一般的に

D = A(T)τ^{n−1}σ′  (2.4)

として表す (Glen の流動則). τ は偏差応力の 2 次の不変量で偏差応力テンソル σ′ を用いて次のように表す.

τ^{2} = \frac{1}{2}tr(σ′^{2})  (2.5)

偏差応力テンソル σ′ と圧力は次のように表す.

σ′ = pI + σ  (2.6)
p = −\frac{trσ}{3}  (2.7)

I は単位行列である.

Msemi-2001-07-11.tex  2001 年 07 月 11 日 (村田泰洋)

式 (2.4) における n の値は応力実験や観測で定められる値であり, 不確定の要素もあるが, 平均値として定数 3 を用いる (Paterson 1994). パラメータ A は温度の関数でアレニウスの式,

A(T) = A_{0} exp (−\frac{Q}{RT} )  (2.8)

を用いる. A_{0} の値は構成する氷の異方性や不純物に依存する. ここでは A_{0} は一定とする. Q はクリープ活性化エネルギー, R は気体定数である.

歪み速度テンソル D を以下のように定義する.

D = \frac{1}{2} (gradv + (gradv)^{T} )  (2.9)

2.3 境界条件

2.3.1 運動学的境界条件

氷体の底面 (z = 0) 運動学的な境界条件 (流体粒子が壁を越えない) は地形的な変化がなく融解も起きないので氷の流出はないので

w(0) = 0  (2.10)

である.

表面 (z = h) では次のように表す.

w(h) = \frac{∂h}{∂t} + v · grad_{H}h + \frac{M_{s}}{ρ}  (2.11)

ここで grad_{H} は

grad_{H} = (\frac{∂}{∂x} , \frac{∂}{∂y} , 0)  (2.12)

である. M_{s} は表面の正味の質量収支である.

2.3.2 力学的境界条件

表面では表面に沿った応力は働かないので表面に対する法線ベクトルを n = n_{i} とすると次のように表す.

(σ_{ij}n_{j} − n_{i}σ_{ij}n_{j})|_{z=h} = 0  (2.13)

\underline{寒冷氷床モデルの定式化  2 寒冷氷床モデルの定式化  5}

氷体全域で圧力融解点以下の温度なので底面では氷体は滑べらず固定されているとする.

u(z = 0) = 0, v(z = 0) = 0 (2.14)

2.3.3 熱学的境界条件底面では地殻熱流量 γ を与える.

k \frac{∂T}{∂z} |_{z=0} = γ  (2.15)

表面温度は関数 T_{s} = T (x, t) として適宜与える.

2.4 直線直行座標での表示

直線直行座標での方程式系は以下のように表現される.

\frac{∂u}{∂x} + \frac{∂v}{∂y} + \frac{∂w}{∂z} = 0  (2.16)
\frac{∂σ_{xx}}{∂x} + \frac{∂σ_{xy}}{∂y} + \frac{∂σ_{xz}}{∂z} = 0  (2.17)
\frac{∂σ_{yx}}{∂x} + \frac{∂σ_{yy}}{∂y} + \frac{∂σ_{yz}}{∂z} = 0  (2.18)
\frac{∂σ_{zx}}{∂x} + \frac{∂σ_{zx}}{∂y} + \frac{∂σ_{zz}}{∂z} = ρg  (2.19)
ρc_{p} \frac{∂T}{∂t} = −ρc_{p}(v · grad)T + k∇^{2}T + ∑_{i,j} D_{ij}σ_{ij}  (2.20)
( \begin{matrix}D_{xx} & D_{xy} & D_{xz} \\ & D_{yy} & D_{yz} \\ & & D_{zz} \end{matrix} ) = A(T)τ^{n−1} ( \begin{matrix}σ′_{xx} & σ′_{xy} & σ′_{xz} \\ & σ′_{yy} & σ′_{yz} \\ & & σ′_{zz} \end{matrix} )  (2.21)
τ^{2} = \frac{1}{2}(σ′_{xx}^{2} + σ′_{yy}^{2} + σ′_{zz}^{2}) + (σ′_{xy}^{2} + σ′_{yz}^{2} + σ′_{zx}^{2})  (2.22)
σ′_{ij} = σ_{ij} + pδ_{ij}  (2.23)
p = −\frac{σ_{xx} + σ_{yy} + σ_{zz}}{3}  (2.24)

Msemi-2001-07-11.tex  2001 年 07 月 11 日 (村田泰洋)

\underline{寒冷氷床モデルの定式化  2 寒冷氷床モデルの定式化  6}
( \begin{matrix}D_{xx} & D_{xy} & D_{xz} \\ & D_{yy} & D_{yz} \\ & & D_{zz} \end{matrix} ) = ( \frac{∂u}{∂x} \frac{1}{2} (\frac{∂u}{∂y} + \frac{∂v}{∂x} ) \frac{1}{2} (\frac{∂u}{∂z} + \frac{∂w}{∂x} ) ; \frac{∂v}{∂y}  \frac{1}{2} (\frac{∂v}{∂z} + \frac{∂w}{∂y} ) ; \frac{∂w}{∂z} ) (2.25)

2.5 浅氷近似

氷床の水平方向と鉛直方向の特徴的長さを比較すると非常に薄い^{1} ので浅氷近似 (shallow ice approximation) を施すことで方程式系は簡略化される (Hutter 1983). 2.5.1 運動量保存則まず, O((H/L)^{2}) の項を無視する近似を施す (Blatter 1995). 式 (2.19) で水平方向の応力の勾配が鉛直方向の勾配に比べて無視できるので

|\frac{∂σ_{zx}}{∂x} | , |\frac{∂σ_{zy}}{∂y} | ≪ |\frac{∂σ_{zz}}{∂z} |  (2.26)

である. 従って, 鉛直方向の運動方程式は次のように表される.

\frac{∂σ_{zz}}{∂z} = ρg  (2.27)

大気圧を 0 とすれば, 表面 (z = h) から z まで積分して,

σ_{zz} = −ρg(h − z)  (2.28)

が得られる. 偏差応力の定義から,

σ′_{xx} + σ′_{yy} + σ′_{zz} = 0  (2.29)

である. 従って, σ_{xx}, σ_{yy} は

σ_{xx} = 2σ′_{xx} + σ′_{yy} + σ_{zz}  (2.30)
σ_{yy} = 2σ′_{yy} + σ′_{xx} + σ_{zz}  (2.31)

となる. 式 (2.30), (2.31) を運動方程式に代入することにより

\frac{∂σ′_{xy}}{∂y} + (2\frac{∂σ′_{xx}}{∂x} + \frac{∂σ′_{yy}}{∂x} ) + \frac{∂σ′_{xz}}{∂z} − ρg \frac{∂h}{∂x} = 0  (2.32)
\frac{∂σ′_{xy}}{∂x} + (2\frac{∂σ′_{yy}}{∂y} + \frac{∂σ′_{xx}}{∂y} ) + \frac{∂σ′_{yz}}{∂z} − ρg \frac{∂h}{∂y} = 0  (2.33)

1実際に存在する氷床の場合, H/L ∼ O(10^{−3})

Msemi-2001-07-11.tex  2001 年 07 月 11 日 (村田泰洋)

\underline{寒冷氷床モデルの定式化  2 寒冷氷床モデルの定式化  7}

が導かれる. 次に O(H/L) の項を無視する近似を行う (Blatter 1995).

|\frac{∂σ′_{ii}}{∂x_{j}} | ≪ |\frac{∂σ_{iz}}{∂z} |  i, j = 1, 2  (2.34)
|\frac{∂σ′_{xy}}{∂x_{i}} | ≪ |\frac{∂σ_{iz}}{∂z} |  i = 1, 2  (2.35)
|σ′_{ii}|, |σ′_{ij}| ≪ |σ_{iz}| i, j = 1, 2  (2.36)

これにより (2.32),(2.33) は

\frac{∂σ′_{xz}}{∂z} = ρg \frac{∂h}{∂x}  (2.37)
\frac{∂σ′_{yz}}{∂z} = ρg \frac{∂h}{∂y}  (2.38)

となる. 2.5.2 エネルギー保存則エネルギー保存則の拡散項および粘性変形項は次のように近似される.

k∇^{2}T ≃ k \frac{∂^{2}T}{∂z^{2}}  (2.39)
∑_{i,j} D_{ij}σ_{ij} ≃ 2(D_{xz}σ′_{xz} + D_{yz}σ′_{yz})  (2.40)

これによりエネルギー方程式は次のように表される.

\frac{∂T}{∂t} = − [u\frac{∂T}{∂x} + v \frac{∂T}{∂y} + w\frac{∂T}{∂z} ] + \frac{k}{ρC_{p}} \frac{∂^{2}T}{∂z^{2}} + 2(D_{xz}σ′_{xz} + D_{yz}σ′_{yz})  (2.41)

2.5.3 構成方程式構成方程式に関係する (2.24), 式 (2.25) は浅氷近似により

τ^{2} ≃ σ′_{xz}^{2} + σ′_{yz}^{2}  (2.42)
D_{xz} ≃ \frac{1}{2} \frac{∂u}{∂z}, D_{yz} ≃ \frac{1}{2} \frac{∂v}{∂z} (2.43)

となる.

Msemi-2001-07-11.tex  2001 年 07 月 11 日 (村田泰洋)

2.5.4 境界条件

力学的境界条件では表面に沿った応力は働かないとした. 浅氷近似の元では式 (2.36) なので表面における応力は

σ′_{xz} = 0,  σ′_{yz} = 0  (2.44)

という条件が与えられる.

参考文献

• Blatter, H. 1995. Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients. Journal of Glaciology 41(138), 333-344

• EISMINT http://www.esf.org/life/lp/EISMINT/lp 008a.htm

• Hutter, K. 1983. Theoretical glaciolog: material science of ice and the mechanics of glaciers and ice sheets. D. Reidel Publishing Co./Tokyo, Terra Scintific Publishing Co.

• Huybrechts,P., T.Payne, and The EISMINT Intercomparison group. 1996. The EISMINT benchmarks for testing ice-sheet models. Annals og Glaciology 23, 1–12

• Paerson, W.S.B. 1994. The Physics of glaciers(3rd ed.). Pergamon, Oxford