夏は暑く冬は寒いことは温帯にくらす人々の常識である。 これは、基本的には、太陽と地球の幾何的位置関係から、 地球のある緯度帯の表面積あたりに到達する太陽放射エネルギーが 年周期で変化していることからきている。ただし、定量的にどれだけの 季節変化があるかを決めるためには、地球から出て行く放射や、緯度帯間の エネルギー輸送もまた季節変化することを考慮しなければならない。
井戸水や湧き水を知っている人ならば、地下水が夏には気温より冷たく、 冬にはあたたかいことに気づいているだろう。温度を年周期で変化させる 原因は地表面から上にあり、地下にその情報が伝わるのは基本的に熱伝導に よっている。地下深くなるごとに、地表面からの影響が弱くなることは たやすく想像できる。しかし、振幅が深さの関数としてどのように減って いくか、また、年周期変化の位相の進み遅れはどうか、といったことに ついては、常識だけでは判断しにくい。そこに、理論に基づく シミュレーションをしてみる価値が出てくる。また、もし実際に測定した 地中の温度分布が熱伝導にもとづくシミュレーションと違っていたとすれば、 シミュレーションにあたって仮定したことのいずれかを疑う必要が出てくる。
地球上の気候の分布の特徴を決める要素として、緯度についで重要なのが、 海陸の分布である。海洋では大陸に比べて気温の季節変化の振幅が小さい。 この原因を簡単に、「海は陸よりも熱容量が大きい」と述べることがある。 これは、「年変化に実質的に効く部分の熱容量が大きい」という意味である。 水が物質として岩石や土壌に比べて比熱(質量あたりの熱容量)が大きい、 ということがその理由ではない。実際、多少大きいが、桁は違わないし、 水のほうが密度が小さいので、体積あたりになおすと違いはますます小さい。 海では水が鉛直方向に動くことができ、物質の移動に伴ってエネルギーが 動くことができる、つまり、熱伝達の種類としての「対流」が起こることが 主要な違いである。対流のおきやすさは、上下の密度差によって変わるので、 熱伝導のように「エネルギー輸送は温度勾配に比例する」わけではないが、 非常に大まかに言えば、熱伝導係数が桁ちがいに大きくなったことに相当する。 その結果、海のほうがより深い所まで年変化の影響がおよぶのである。
海と陸の気温の年変化の振幅の違いはよく知られているが、位相の違いは どうだろうか。
気象観測点の、毎日の気温データに、 最小2乗法で1年周期のサインカーブをあてはめてみた。 ([気温の年変化の振幅と位相]参照)。 タクラマカン砂漠と日本では、気温が極大になる時期が約1か月違う。 日本文化の中には中国からきた暦が生き続けているが、「立春」「立秋」は それぞれ寒いさかり、暑いさかりであり、名前が実感に合わない。 しかし、唐の長安(現在の西安=シーアン)では、確かに気温の極小・極大を 過ぎて、春や秋への動きを感じさせる時期にちがいない。 海に近い地点の気温には、海上からの移流の影響があることを考えると、 これは海と陸との「年変化に効く熱容量」の違いで説明できる。
陸または海の表層のエネルギー保存を考える際には、 有限の厚さをもった層について考える場合と、 「地表面」あるいは厚さを無限に小さくした層について考える場合がある。 [気候システム論:陸と海の地表面のエネルギー収支の季節変化] [HARTMANN 4章(学内限定)] ここでは、有限の層を考えてみる。 実際には、地温の年変化の振幅は地面から下に向かってだんだん小さくなっていく。 それを熱伝導の式に基づいて扱うのは次回の問題とし、 ここでは、陸では1メートル程度の厚さの層が一様な温度をもつという 板のような近似で考えてみよう。 海では、鉛直には熱伝導だけでなく対流によるエネルギー輸送があるので、 (対流の強さは実際には温度や塩分の成層状態によって違うが) 100メートル程度の厚さの「混合層」が一様な温度をもつという近似ができる。
土壌または水の層が、温度Tに比例するエネルギーC Tをもつとする。 (その比例定数Cを熱容量という。) 水平方向には同じものがずっと広がっているように考えるのだが、 便宜上、単位面積あたりで考え、Cは単位面積あたりの熱容量とする。 式は次のようにたてることができる。
C dT/dt = Gtop - Gbottom ... (1)
Tは温度で変数、Cは面積あたりの熱容量で定数とする。 fは 時間あたり・面積あたりのエネルギーの出入り(エネルギーフラックス密度)である。 Gtop、Gbottomはそれぞれ層の上側と下側の面を通る 正味下向きエネルギーフラックス密度とする。 陸では土壌層の熱伝導だけを考えればよい。海などの水体では対流もある。 ここでは、板を適当に厚くとればGbottomは無視できると仮定し、 GtopをGと書くことにしよう。 地表面という面自体は(無限に薄い層と考えてもよいが)エネルギーをたくわえることが できないので、Gは、地表面のすぐ上での正味下向きエネルギーフラックス密度に 等しくなければならない。後者をやや具体的に表現すれば
G = Rnet - H - L E ... (2)
ただしRnetは地表面での正味下向き放射、 Hは正味上向き顕熱輸送、 L Eは正味上向き潜熱輸送(Eが蒸発の質量輸送、Lは水の質量あたりの蒸発の潜熱) である。 なおRnetは、太陽(solar)放射・地球(terrestrial)放射の下向き・上向きの4成分に分けられる。
Rnet = RS↓ - RS↑ + RT↓ - RT↑ ... (3)
現実には、(2)の右辺を構成する各項、 たとえばH、L E、- RT↑は、 地面温度と連立された形で決まる(原因と結果を分けられない)。
しかし、まず、 (2)の右辺(Gに等しいが仮にfと呼ぶことにする)を モデルの外部から与えられた強制項(forcing、力ではないが「外力」ともいう)とし、 時間の関数として正弦波(サインカーブ)的に変化するものとする。 具体的には(Fortranの代入文の形で書けば)
f = famp * cos ( 2 * pi * ( t - tphase ) / tcycle ) ... (4)
のようにし、tcycleに周期、fampに振幅(最大と最小の差の半分)、 tphaseには位相を、uが極大となる時刻という形で与えることにする。 piは円周率で、あらかじめ値を入れておく (たとえば pi = 4.0d0 * atan(1.0d0); atanはtanの逆関数)。 プログラム例をFortranとC言語で用意した。
| Unix用(改行コードLFのみ) | Windows用(改行コードCRLF) | |
|---|---|---|
| Fortran | cycle.f | cycle.f |
| C | cycle.c | cycle.c |
コンパイルして実行を開始すると、 端末画面に英語(の断片)で質問が出るので、 それに答える形で、熱容量と、強制項fの振幅の数値を入れると 実行が始まって標準出力に結果が出る。 結果をファイルに入れるには実行開始の際に出力redirection機能(>)を使って ファイル名を指定しておく。 ここで、量の単位は、国際(SI)単位系を使うことを想定した。 時間の単位は秒である。
差分計算の時間間隔やシミュレーション対象の時間の長さなどは、 前の実習課題1では 実行時の入力データとして与えるようにしたが、 今回は、ソースプログラム中に書きこんであるので、 変更したいときは、ソースプログラムを編集してコンパイルしなおす必要がある。
面積あたりの熱容量Cを大まかに見積もってみる。
これを実行してみると、 熱容量と強制項の振幅を変えたことによって 温度の振幅がどう変わってくるかについては結果はもっともらしいのだが、 温度の位相に関してはおもしろくない結果になる。 そこで、モデルを練りなおすことにする。
(4)式のfをf = f1 + f2のように分け、
f1のほうだけを外部からの強制効果と考えることにする。
これは、正味下向き太陽放射の年平均値からの差に対応すると考える。
f1 = famp * cos ( 2 * pi * ( t - tphase ) / tcycle ) ... (5)
太陽放射の季節変化を1年周期の正弦波で近似することは、 温帯ではもっともである。 地表面に達する太陽放射は大気中での吸収・反射によっても変わるが、 大気上端に入射する太陽放射は幾何学的計算で求めることができる。 これを周期解析してみると、中緯度の場合は、分散の90%以上が 1年周期成分に含まれている。 低緯度と高緯度の場合は半年周期成分も無視できない。 ([気温の年変化の振幅と位相] の中の [日射量の年周期成分]参照)。)
f2で表現されるのは、 正味下向き地球放射、蒸発に伴う潜熱輸送、顕熱輸送である。 いずれも、必ずしも比例関係ではないが、 地表面温度が高いほど逃げていくエネルギーが大きい傾向がある。 その特徴を簡単な形で表わすために、 「出ていくエネルギーは、現在の温度と、あらかじめ与えられた一定の温度 (「参照温度, reference temperature」と呼んでおく)との差に比例するとする。 このように、時間変化が一定値との差に比例するような作用を、 「緩和, relaxation」と呼ぶことがある。 また、特に、エネルギーの出入りに関しては、 「Newton冷却, Newtonian cooling」と呼ぶこともある。 ただし、現在の温度が参照温度よりも低い場合も、比例定数は同じだから、 そのような場合は、実際は冷却ではなく加熱の働きをする。 プログラム上の形は、 参照温度をtmpref、比例係数をccoolとすれば、 次のようになる。
f2 = -ccool * (temp - tmpref) ... (6)
プログラム例をFortranとC言語で用意した。 便宜上、tmprefは255 [K]で固定し、 ccoolを標準入力(「<」でredirectしなければキーボード)から入力する データに追加した。
| Unix用(改行コードLFのみ) | Windows用(改行コードCRLF) | |
|---|---|---|
| Fortran | cycle2.f | cycle2.f |
| C | cycle2.c | cycle2.c |
このモデルで時間発展を長く計算すると、温度は1年周期でほとんど同じ値が くりかえすようになる。 その平均値は、あらかじめ与えた参照温度tmprefに等しくなる。 現実の気候とある程度対応がつくのは、熱容量C (プログラムではcapac) を変えたときの、温度の1年周期変化の振幅と位相である。
同じ問題を、数式変形による理論計算で考えてみる。 前の節の設定を含めて方程式(1)を修正すると次の形になる。
C dT/dt = F1(t) - r × ( T - Tref ) ... (7)
現実にはエネルギーの水平移流があるが、 ここでは緩和にむりやり含めていることになる。
強制項は1年周期の正弦波 (t = 0 で F1が極大)と仮定する。
F1(t) = AF cos (ω t ) ... (8)
ω = 2 π/tyear ... (9)
(ただしtyearは1年を表わす定数)。
この方程式は「線形」なので、 強制項の形がsin, cosであるとき、結果も同じ周期のsin, cosだけからなると 考えてよい。 そのような形を想定して方程式(7)に入れ、 係数を求めて整理してやると、解は次の形に書ける。
T(t) = Tref + AT cos ( ω t - φ ) ... (10)
位相: φ = arctan ( ω C / r ) ... (11)
振幅: AT = (AF / r) cos2φ ... (12)
次の図には、 熱容量Cを横軸にとり、緩和の係数rの値をいくつか変えて、 振幅ATを示した。 ただし、便宜上AF = 100 W / m2としてある。 AFを変えた場合、ATはAFに比例して変わる。 Cの小さい極限では、AT = AF/ r。 また、このAFのもとで、 ATが現実にありそうな値 10 K になるような rの値は、10 W m-2 K-1である。
また、次の図には、 熱容量Cを横軸にとり、緩和の係数rの値をいくつか変えて、 温度の1年周期の位相を、Fの極大からの日数に換算した形で示した。 Cの小さい極限では位相差は0、 Cの大きい極限では位相差は90度(91.5日)となる。 振幅が現実と対応するようなrの値を使うと、 位相はCが106 -- 109 SI単位の範囲で 熱容量に敏感である。
さて、式(11)と(12)は、 ATとφを求める形をしているが、 変形すれば、 AF, AT, φの情報からC, rを求めることができる。
C = AF/(2 AT ω) sin 2φ ... (13)
r = AF/(2 AT){ cos 2φ + 1 } ... (14)
細かいことを言うと、この節の問題は、 前の節の時間発展問題と厳密に同じではない。 ここでの解きかたは、 無限の過去から同じ状態が周期的にくりかえしているという 仮定が含まれている。 一方、前の小節で紹介したプログラムは、この仮定とは違う 勝手な初期条件(時刻t = 0で温度T = 288.15 K = Tref + 33.15 K) から出発していた。 しかし、この方程式の場合、勝手な初期条件から出発しても、 長い時間続ければ、ほぼ周期的な変化に落ち着き、 それは初めから周期変化を仮定して得られた解と、事実上同じである。
なお、上の2つの図で使った数値は、 Fortranプログラムamppha.fで計算した。 (これはこの図のための臨時のもので、汎用性は考えていない。 変数の宣言を省略してしまっているなど、模範的でないプログラムである。) 作図の準備として、結果を「amppha.out」というファイルにしておく。 作図は、gnuplotで、 graamp.gpltと grapha.gpltをloadして、PostScriptファイルを 作り、 ImageMagickのconvertコマンドで画像ファイルに変換した。 (この段落のリンク先のファイルはLinux上の作業用のものなので 改行コードはLFのみになっている。 Windows上でうまく読めない場合は個別に相談してほしい。)
今回の演習問題では仮想的に強制項を想定して考えたが、 もう少し現実に即した扱いをするための考えかたを示しておく。
Rnet = H + L E + G ... (2')
R↓は次のように
太陽放射と下向き地球放射(大気からの放射)を合わせたものである。R↓ = RS↓ - RS↑ + RT↓
R↓ = σTs4 + H + L E + G... (2'')
R↓effとしておく)は、
正味放射量の上向き地球放射の代わりに(地表近くの)気温Taの黒体放射を含めたものである。R↓eff = R↓ - σTa4
R↓eff = σ×(Ts4 - Ta4) + H + L E + G... (2''')
Ts - Taがあまり大きくないとき、
σ×(Ts4 - Ta4)は
4σ Ta3×(Ts - Ta)と近似できる。
また、バルク伝達式を使えば、HはTs - Taに比例し、
L Eはqsat(Ts) - qaに比例する(qは比湿)。
地中熱伝導GはTsともっと深い層の地温の差に比例すると考えられる。
つまり(2''')の右辺の項はいずれもTsを緩和するように働く。http://www.asahi-net.or.jp/~rk7j-kndu/に
正誤表がある。]http://www.asahi-net.or.jp/~rk7j-kndu/に
正誤表がある。]
[2009年4月補足:
付録のBasic言語によるプログラム例をFortranで書きなおしてみたので、
ほしいかたは増田までお問い会わせください。]