環境システム科学演習− MATLABを用いた海洋情報処理・作図演習


目次    授業ホーム    ホーム


3) MATLABを用いた海洋観測データ処理(海水物性)・作図

MATLABは多くの海洋学者が利用している為、海洋学専門のプログラムがインターネット上で無償で配られている。

http://woodshole.er.usgs.gov/operations/sea-mat/

もその一つである。


-1 塩分


目次    授業ホーム    ホーム

海水には、様々な物質が溶け込んでいる。その主なものは、Chlorine ion (塩素イオン)Sulphate ion (硫酸イオン)Sodium ion(ナトリウムイオン)、Magnesium ion (マグネシウムイオン)Potassium ion (カリウムイオン)等である。海水にはこれらの溶存塩総量が場所によって様々な濃度で溶けているが、興味深いことに世界中の海で、ほぼ一定の相対的な割合で溶け込んでいる。これは、長い地球の歴史の中で如何に海水がよく混ぜられているかを表している。1kgの海水中に、何gの溶存物質が溶け込んでいるかで定義されるのが塩分である。具体的には、全ての海水中の炭酸塩を酸化し、ヨウ素と臭素を塩素で置き換え、全ての有機物質を酸化した際の固形物質が海水1kg中に何gであるかで定義された。この塩分は外洋域では大体35ppt(parts per thousand)程度となる。上記の方法は、繰り返し分析するのに非常に面倒であった為、Clorinity(塩素量)のみを硝酸銀滴定により測定し、塩素量と溶存塩総量との経験式(塩分=1.80655×塩素量)から塩分を求める方法を便宜的にとった。

 しかしながら、現在では塩分は電気伝導度と水温、圧力の観測を同時に行い、経験的に求める方法をとる。これは、電気伝導度を用いた方が、塩素量を分析するよりも実用的に精度を一桁向上させるためである。このため、電気伝導度でもとめた塩分を(practical salinity)実用塩分と呼び、塩分35pptと同等な塩分は塩分35と表記するか、あるいは極稀に塩分35 psu (practical salinity units)と表記される。この実用塩分は水温15℃、1気圧下では以下のような経験式として表される。

ここでSは実用塩分、K15は電気伝導度率とよばれ、塩分Sで水温15℃、水圧0の時の電気伝導度C [mS/cm]の、塩分35のときのそれに対する比として定義される。即ち、K15=C(S, 15, 0)/C(35, 15, 0)。定数係数は、

a0

0.0080

a1

-0.1692

a2

25.3851

a3

14.0941

a4

-7.0261

a5

2.7081

である。しかしながら上記では、CTDなどの測器で連続的に測定される異なる圧力や温度下における電気伝導度と塩分の関係を記述することができない。そこで、水温と圧力の電気伝導度に及ぼす影響を考慮できる経験式が導かれるに至った。まず、測定した電気伝導度C(S,T,P)を塩分35水温15℃、大気圧の電気伝導度で比を取りRとする。即ち、R=C(S,T,P)/C(35,15,0)。この比は以下のように分解される。

このうちRTは大気圧下で任意の水温、塩分の海水の電気伝導の大気圧下、塩分35の海水のそれとの比である。観測した任意の水温、塩分、圧力下の電気伝導度を(6)式を用いてRTに変換し(RTとは定義から、RRprTで割ったものである)、RTについて経験式をつくって水温と圧力の効果を考慮する。過去の室内実験からRprTの経験式が以下のように求められている。

A1

2.070 x 10-5

A2

-6.370 x 10-10

A3

3.989 x 10-15

B1

3.426 x 10-2

B2

4.464 x 10-4

B3

4.215 x 10-1

B4

-3.107 x 10-3

 

c0

6.766097 x 10-1

c1

2.00564 x 10-2

c2

1.104259 x 10-4

c3

-6.9698 x 10-7

c4

1.0031 x 10-9

求めたRTを用いた経験関数として塩分は以下のように近似できる。これらは過去の大気圧下の室内実験から求められたことから、RTに変換する必要があるわけだ。T=15℃とすると、(9)式は(4)式と同一である。

b0

0.0005

b1

-0.0056

b2

-0.0066

b3

-0.0375

b4

0.0636

b5

-0.0144

k

0.0162

C(35,15,0)

42.914

塩分を現場水温、電気伝導度、圧力の関数として求めるMATLAB関数プログラムは以下のようになる。

上記プログラムは、functionから始まるMatlab M-Fileで、関数として機能する。関数内で計算に用いられる変数は、関数内のみで有効であり、function行において、等号=の左辺で定義した以外の変数はワークスペースに出力されることはない。このように決まった計算を繰り返し行う必要があるようなケース(例えば塩分の計算)では、計算式をあらかじめ関数プログラムとして記述しておく。これによって、何度も行う作業部分を他のプログラムから簡単に呼び出して使うことが出来る。

ピリオドを3連続で打つと(...1行を改行してプログラムを書ける。改行後と前とはプログラム的にはつながっている。

 


課題 塩分の計算と図示


目次    授業ホーム    ホーム

以下のリンクから電気伝導度と水深データをダウンロードしよう。電気伝導度データを読み込み、塩分に変換して塩分の鉛直分布図を書いてみよう。

TempCondPress.mat


-2 水温、現場水温 ポテンシャル水温


目次    授業ホーム    ホーム

水温は圧力によっても変化する。即ち、圧力変化によって海水がほんの少しだけ圧縮されたり膨張したりする際に、気体と同様に加熱、冷却する。深海での高圧力のため、深海における水温が若干上昇し、海水の縦方向の安定性を考える際に誤解を生じやすい。

図3−1

このため、水温の定義には現場水温と、ポテンシャル水温が存在する。現場水温は、たとえば、水温センサーを降下しその深さで計測された水温のことである。ポテンシャル水温は、海水を断熱的に(まわりとの熱のやり取りを遮断し)海面まで上昇させた時の海水の温度を指す。ポテンシャル水温を計算することによって、圧力による温度上昇分を取り除いて議論できる。上図は圧力によって現場水温が上昇している例である。一方ポテンシャル水温Θは下向きに上昇していない。

今、水圧poにある塩分So、水温toの海水を圧力prまで断熱的に上昇、下降させるときの温度を考えると以下のように書ける。

ここでpr0とすると、Θはポテンシャル水温である。右辺第2項では、単位圧力当りの水温変化率Γを圧力についてpoからprまで積分し、その分水温の初期条件toから変化が生じることとなる。このΓは、圧力をpoからprに連続的に変化させ水温が断熱的に変遷する中で変化する。この積分記号内の関数Γも結果として得たいΘの関数である為、数値的な積分方法によってはdpを大きくとった場合無視できない誤差を生じてしまう。コンピューターが現在の性能に満たなかった頃には、Γを効率的に数値積分する方法が考案された。数値積分には4次のRunge Kutta法が用いられた(Fofonoff 1977)。下記のプログラムは、4次のRunge Kutta法を用いたものである。Runge Kuttaについてはこの授業でカバーしないが、興味がある人は、このリンクに説明をつけたので参考までに。このプログラムを用いて、塩分の計算に用いたファイル中の水温をポテンシャル水温に変換しよう。ここで、theta(s,t,p,pr)のs,t,pは現場の塩分、水温、圧力(dbar)である。prは基準面で、pr=0にとれば計測した水を断熱的に海面に上昇させたときのポテンシャル水温がptとして出力されることになる。プログラム中で使用している関数atgは、上記で述べたΓである。Γは塩分、水温、圧力の関数である。Γには経験式が存在し、function atgとしてsample5に示してある。これらのサンプルプログラム4、5をfunction名で保存しよう。

 


課題 ポテンシャル水温の計算と図示


目次    授業ホーム    ホーム

上記プログラム作成後、これらのfunctionを用いてTempSaltPressAtlantic.mat中の現場水温をポテンシャル水温に直し、現場水温の鉛直分布と比較する図を書いてみよう。

 


-3 密度 ポテンシャル密度 σ


目次    授業ホーム    ホーム

水は大氣と違いほとんど圧縮できない流体である。例えば注射器に水を満たし、先端をしっかり止め、ピストンを押しても全くと言ってよいほど、ピストンを押して水を圧縮することは出来ない。しかしながら、深海における水圧は大変大きい為、しばしば圧力による圧縮効果を無視することが出来なくなる。例えば、塩分35、水温0℃で水面にある水の密度が、1028.13 kg/m3であるのに対し同じ塩分,水温で水深4000mにある場合だと、密度は1048.49 kg/m3である。圧力によって縮められる分、密度が2%程度大きいことになる。この圧力による密度の変化、即ち水の圧縮性について圧縮率βを以下のように定義する。

βの逆数としてKは以下のように書いても良い。

室内実験でKを求める為に(12)式の-Vdp/dvで表される勾配を、圧力を0からp'まで変化させた時の平均圧力pと、はじめは体積がVoであったものが圧力によってVに変化したとして、簡略化して以下のように表した。

圧力を受けた前後で、海水の質量自体は変化しないことから、Vをその質量で割って海水1kgがしめる体積、比容に直して記述しても差し支えない。このKをpについての2次関数、Ko(T,S)を水温、塩分の関数として室内実験から求めた。A,Bは水温の変化に伴う効果を含めて求められた。

Kに関する経験式は以下のように求まった。

(14)式から、圧力によって変化する比容、密度を記述する式を与える。

(16,17)式中のαo、ρoは、それぞれ大気圧下の海水の比容および密度である。海水の重さは、圧力以前に現場の水温と塩分の関数であることは言うまでもない。暖かい(冷たい)海水ほど軽く(重く)、塩分が高い(低い)海水ほど重い(軽い)。圧力を大気圧に設定し、水温と塩分の変化に対する海水の密度変化を調べる室内実験も行われ、その経験式は以下のように記述される。

さて、海水の密度について述べてきたが、海水の密度の変化は大きくても数%程度である。前述の例だと、海面で1028.13 kg/m3であった海水は、水深4000mの圧力のみの圧縮効果で1048.4913 kg/m3に増加する。これはわずか全体の約2%の増加にすぎない。このように海水の密度は大きく変化したとしてもせいぜい下二桁で変化が起こる。そこで密度から1000を単純に差し引いて変化を判り易く書き記すのが慣例となっており、密度から1000引いた量をσと定義している。

σはもちろん水温、塩分、圧力の関数σ(T,S,P)である。σには様々な種類が有る。

 

まず上記のσは現場水温、塩分、圧力の関数で現場の海水密度を表す。

σTは、σ(T,S,0)と圧力を大気圧と仮定し現場の水温、塩分から(18)および(19)を用いて求められる。σTは水平的な海水の密度変化を示したいときに頻繁に使われている。


-4 ポテンシャル密度 σΘ


目次    授業ホーム    ホーム

上記のσTは圧力の効果を考慮できないため、比較的深い水深で用いると、そこでの圧力で上昇した現場水温の影響を受けることになる。これは、もともと圧力を考えないσTの概念からして想定外となる。またσは圧力の効果を加味し現場の密度を知ることが出来るものの、前述のとおり圧縮自体の効果と圧縮による水温上昇効果が含まれており、圧力が増減する縦方向に密度を比べることは好まれない。そこで、圧力にともなう水温増減効果、圧縮による加密化を取り除いて縦方向の水塊比較を行う為に、前述したポテンシャル水温を用いて水面における密度を計算することが従来行われている。即ち、海水を断熱的に海面に上昇させ、圧力による水温上昇効果と、圧縮自体の効果を取り除いた上で密度を計算する。これは、異なる水深に位置する異なる海水を同じ基準面に並べ同じ水圧下ならば2つの水塊がどう比較できるか確かめることにあたる。ポテンシャル密度の場合、基準面は水面である。この方法で得られる密度をポテンシャル密度と呼ぶ。またρθから1000引いた値をσθと呼ぶ。即ち、

このσθをσTと比較した例は、図3−1に示してある。σTを現場水温から算出した場合、密度が下層で小さくなる。即ち、下層の方が軽く不安定であるかのように見える。しかし、σθをポテンシャル密度から計算した例では下層は上昇よりも重たいことがわかる。

しかし、σθも以下の図3−2に示す太西洋深海の例では万能ではないことがわかる。すなわち、海底付近のσθは下層のほうが小さく不安定となってしまう。この密度逆転構造は当初議論されたが、実はポテンシャル密度の計算があまりにも深い水深の海水には妥当ではないことによる。実際はここでの塩分は海底より若干上で高く現場水温が高い。その下層では塩分が低く水温が低くなっている。密度は少しだけ下層の方が大きい。ところで圧縮率は海水温によって変化する。冷たい海水の方がより圧縮性が高いことから、下層の低温低塩水は断熱的に5000m上の海面まで持ち上げるとより膨張し、ポテンシャル密度がより下がる。それより上の高温高塩水は、下層と比べて高温であるため膨張が比較的盛んでなく、ポテンシャル密度はさほど小さくなれない。このため結果として、ポテンシャル密度の逆転が起こっている。

図3−2

問題は、海面まで断熱的に持ち上げる際の膨張効果が大きすぎる為である。例えば、基準面を海面ではなく水深4000mに取った場合この密度の逆転は解消されることがわかる(図3−3)。

図3−3


課題 ポテンシャル密度(σθ)の計算と図示


目次    授業ホーム    ホーム

下記リンクから、海水のσ(T,S,P) (現場密度-1000)を求めるプログラムをダウンロードし、

1.それを用いてTempSaltPressAtlantic.mat中のσTおよびσを求めるプログラムを作成し、

2.両者を重ねて図示し比較しよう。

次にsigma_p.m等を用いてσθを求めるプログラムを作成し作図もしよう。このとき水深4000-6000mの分布に着目しよう。

4.また基準面水圧を水面0でなく、4000dbarとしσ4を求め図示しよう(このときも深4000-6000mの分布に着目しよう)。

σTを求めるプログラムは(18)式から自分でfunctionを作ってもよいし下記のリンクのものを使っても良い。

σ(T,S,P)をもとめるプログラム