﹂ エ

淳 エ

﹁

jL

方

﹁

d

一 副

散ム 仏

﹂1 守 d

流移

ここで，u ( X， y， t  ，)

V  ( 

X， y， t )は x，y方向の水深平均流速，

c  ( 

^{X， y， t )は潮位，}

c  ( 

X， y， t)  は拡散物質の水深平均濃度，q'( X， Y )は流入水量，

s  ( 

^X， Y )はソースもしくはシンク等 を表す項，/はコリオリ係数，九 2は海底摩擦係数(以下， 0.0026と仮定した) ， V

，  ( 

^X， 

y)は水平方向の渦動粘性係数である.

計算の手順は，渦動粘性係数 vIが 未 知 量 で あ る た め ， 先 ず ¥ を 一 定 値 と し て 適 当 な 値を仮定した.次に，

( 4 . 3 )  ‑ ‑( 4 . 5 )

式を用い有限差分法により

M

₂潮(周期 :T= 12時間 25分 ) を 対 象 潮 汐 と し て 潮 流 計 算 を 行 い ， 場 所 毎 の 潮 流 最 大 流 速 に を 求 め る .

( 4 . 3 )   ‑ ‑

(4.5)式の差分化は日野・伊藤(1971)を参照し，慣性項に l次精度の風上差分を，粘性項 に

2

次精度の中央差分を採用した.ここで 渦 動 粘 性 係 数 が 分 散 係 数 D と近似的に等し いと仮定し，得られた V^Iと (4.2)式 よ り 場 所 毎 の ¥ を 推 定 す る . 新 し く 得 ら れ た 円 を 用いて再び潮流計算を行い 新たに得られた Vmに よ り ¥ を 再 度 推 定 す る . こ の プ ロ セ ス を 繰 り 返 し 行 う こ と に よ り ¥ の 収 束 値 を 求 め た . 最 終 的 に 得 ら れ た ¥ の 収 束 値 を 用 いて潮流計算を行い その

1

周期分の計算結果を水質拡散シミュレーションに用いた.

水質拡散シミュレーションでは，保存拡散物質として塩素イオン濃度(東京湾におい ては塩分濃度)を対象とし，場所毎の分散係数 D には最終的に得られた

v t

の値を用い て拡散計算を行った.塩素イオン濃度や塩分濃度は，内湾においては河川などから流入 した淡水による希釈のため一般的には湾口から湾奥にかけて減少していく濃度分布を持 つ.計算方法は，保存形である (4.6)式を対流形に書き直し，次式のように splitoperater  approach  (もしくは，時間分割法 (time‑splittingmethod)  [Yanenko (1971)参照]とも呼 ぶ)を適用しx，y方向の移流と拡散を分離して計算を行った.

dC  . 

^{T T}  

dC 

~+U 一一 =0

d t   d X  

^(4.7a) 

(4.7b) 

子市(ま [(h+C)Df14(h+5)Dt+s‑q c ) μ  

移流項の計算には Komatsuet aJ. (1985)により開発された特性曲線法に基づいた 6‑point schemeを ， 拡 散 項 に は 有 限 差 分 法 (2次精度中央差分)を用いることにより一周期の平 均濃度が一定に収束するまで計算を行った.なお， 6‑point  scheme は Falconerand Liu  (1988) に よ り 行 わ れ た 数 値 実 験 に よ り ， 一 般 的 に 精 度 が 良 い と 言 わ れ て い る QUICK schemeより更に精度が良いことが確認されている.

各湾においてこの一連の計算を種々の比例定数

α

について行い，その計算結果と実測 データを比較検討し 最も現況をよく再現している

α

をその湾における比例定数とし た 潮 流の 計 算 結 果 は 渦 動 粘 性 係 数 の 変 化に 対 し て 応 答 が 鈍 くほとんど変化しないこ と，また潮流の実測結果に高い精度が期待できないことから，比例定数

α

についての検 討には拡散シミュレーションの計算値と塩素イオン濃度(もしくは塩分濃度)の実測値 の比較により行った.拡散シミュレーションの計算結果は分散係数の変化に敏感に応答 するため， αの決定に際し高い精度が期待できる.

4.2.3計算結果について

各湾における計算条件を表‑4.1 に示す.潮流シミュレーションにおける湾口での境 界条件は， M2j朝を対象潮汐にして周期 T=12時間 25分，振幅 αの正弦波を与えた.各 湾の湾口での潮位振幅αの値は潮位表[気象庁(1988)Jより一年間の干満差のデータか ら年平均振幅を算出し適用した.コリオ リ係数は湾内で一定とする f‑plane近似を採用し た.各湾の地形データ(平均水深 ・海岸線の地形)は海図より読み取り，淡水流入量の データは各湾の管理者の調査結果[博多湾:博多港港湾管理者(1989)，伊勢湾:三重県 (1988)，東京湾:千葉県(1992)・神奈川県 (1992a)， 大 村 湾 : 長 崎 県 衛 生 公 害 研 究 所 (1992)J より，またその他不明な河川流量については流量年表[日本河川協会 (1989)J

より一級河川の流量と流域面積から推定した値を採用した.潮流シミュレーションの計 算結果は， 実測結果[博多湾:博多港港湾管理者(1989)，有明海 :海上保安庁 (1974)， 鹿児島湾 :海上保安庁(1979)，伊勢湾:日本海洋学会沿岸海洋研究部会(1985b)，東京 湾:同左(1985a)，大村湾:同左(1985c)]と比較し再現性の確認を行った.

塩素イオン濃度の拡散シミュレーションでは，湾口部における境界条件として潮流が 湾内に流入するときには実測結果[博多湾:博多港港湾管理者(1989)，有明海 :気象庁 (1974)，鹿児島湾:鹿児島県(1989)，伊勢湾:三重県(1992)・愛知県(1992)，東京湾:神 奈川県(1992b)・東京都(1991)，大村湾:長崎県衛生公害研究所(1992)]をもとに設定 された一定濃度を与え，流出時には自由流出条件 (dC/dn 

= 

0 : nは湾口境界に直交する 方向)を与えて計算を行った.また，計算結果との比較にも同じ実測結果の年平均値を 用いた.

博多湾における潮流シミュレーションの計算格子を図‑4.2に 最大上げ潮時と潮汐 残差流の計算結果を図‑4.3a)，b)に，塩素イオン濃度拡散シミュレーションの結果を図

‑4

.4に示す.また，計算結果と実測値の比較を容易にするために 湾口から湾奥にか けての一次元的な濃度分布により計算結果の比較を行った.その結果を図‑4.5 に示 す.比較地点は図‑4.4中に黒点で示している通り，湾口から湾奥にかけての 7点であ る.図‑4.5より博多湾における比例定数の最適値は α =90程度であることが分かる.

α= 90のときの分散係数の分布を図‑4.6に示す.

同様に有明海，鹿児島湾，伊勢湾 東京湾 大村湾の

5

つの湾についても拡散シミュ レーションの計算結果と実測値との比較検討を行った.その結果 各湾における比例定 数

α

の最適値は，有明海で α=15程 度 鹿 児 島 湾 で α=35程 度 伊 勢 湾 で α=330程 度，東京湾で α=50程度，大村湾で α=200程度をとることが分かった.有明海，鹿児 島湾についての計算結果と実測値との一次元的濃度分布の比較を図‑4.7，4.8に示す.

また，伊勢湾の塩素イオン濃度分布の計算結果を図‑4.9に 東京湾の分散係数分布の 計算結果を図‑4.10に，大村湾における塩素イオン濃度分布の計算結果を図‑4.11 に それぞれ示す.

4 . 3 比例定数 α の普遍関数表示化

各湾において最適値として求められた比例定数

α

は それぞれ異なる値をとった

.α

の値に対する影響因子は様々なものが考えられるが 各湾における成層度や吹送流によ る影響が最も強く分散効果に寄与しているものと考えられる.そこで 成層度と吹送流 による効果が湾全体に与える影響を表す overallなパラメータを用いて

α

の普遍関数表 示を試みる.

4.3.1成層度パラメータに関する考察

対象水域全体の成層化の度合を表すパラメータは従来より種々考えられている.例え ば，Harleman (1966)による stratificationparameterや Harlemanand Ippenによる estuary

number  [Fischer (1972)参照]， Fischer (1972)による estuarineRichardson number，また 内部 Froude数等もあげられる.室内実験結果をもとにした Fischer(1972)による解析に よると非成層場における分散係数 D。と成層している場の分散係数 Dの比 D/D。は内部 Froude数にはほとんど依存せず ，estuarine Richardson numberに強く依存していることが 報告されている.また， stratification parameterと estuarynumberは estuarineRichardson  numberにより表現できることも記されている.

これらのことより，本研究では成層度を表すパラメータとして estuarineRichardson  numberを採用した.湾全体の成層度パラメータとしての estuarineRichardson number' 'P 

は次式の様に定義される.

ψ ‑ A   (4.8) 

ここで，.1

p  /  P  ( 

= 0.02 )は淡水と海水の相対密度差，Qfは単位時間当たりの淡水流入 量，Bは湾口の幅，V

Tは湾口における潮流流速の l周期間の rootmean squareである.

α

とVの評価に必要な各湾における諸量，及び Vの値を表‑4.2に示す.また，

α

と Vとの関係を図‑4.12に示す.次に述べる吹送流パラメータが各湾で一定でないため，

α

にばらつきはあるものの，ほぼデの 1/2乗に比例していることが分かる.

4.3.2吹送流パラメータに関する考察 ^ノ

風により生じる吹送流が分散能に与える影響を表すパラメータ φ を次式の様に定義

した.

φ _pG 戸 U~

^Abay 

p  Qr  Vf 

^(4.9) 

ここで，

P n  

= 1.293 ( kg/m³ )は空気の密度，

P

は海水の密度，

u w

^{は湾における年平均風}

速，

y=

而 茄26は水表面上での風の摩擦速度 UuとUwの比，A_bay は湾の表面積 ，^{Q T}は潮 汐により湾口から湾内へ流入する入潮量の半周期平均値である.φ は風によって生じる 水表面でのせん断応力が湾全体に及ぼす仕事量

pafuf‑γuw‑Ab

、と潮汐が湾に与える 運動エネルギ‑ pQT VT2の比として表現されている.

φ

の評価に必要な諸量及び φの値を表‑4.2に示す

.uw

については各湾周辺の気象台 などの観測結果より各湾の代表的な地点における年平均値を採用した

.α

とφの関係を 図‑4.13に示す.デが一定でないため

α

はばらつきを示しているが，

α

は見かけ上 φ のほぼ 1/3乗に比例することが分かる.

『剛"，.‑

4.3.3 tpとφによる比例定数

α

の関数表示

比例定数

α

が 個 別 に は そ れ ぞ れ ヂ1/2とφ1/3にほぼ比例していたことから，

α

とそれ らの積 tp1/2・φlβ との関係を求めることを試みた.

試行錯誤の結果より比例定数

α

は，ヂ1/2・φ1/3のほぼ 1/2乗に比例することが分かっ た.したがって，分散係数Dは Vとφにより次式のように定式化されることになる.

D  =  7 0   p

1l4

φ

1/6 

V 

m 

h 

(4.10) 

図‑4.14に

α

とtp1/4・φ1/6の関係を示す.図 ‑4.14は，図 ‑4.12，4.13に比べてばら つきが小さくなっており，このことから淡水流入による成層化と吹送流が年平均的な分 散能に最も支配的な影響を与えているものと推測される.

この結果より，式 (4.8)，(4.9)， (4.10)を用いることにより平面 2次元モデルのための分 散係数を精度良く推定できるものと思われる.

4 .4分散係数の異方性についての検討

平面

2

次元モデルにおける移流拡散方程式 (4.6)は，一般的には次式の様に書ける.

三 d t  ¥  (HC) ‑ ‑

^I 

+ 立川町)+立 d x  ¥ 

^I 

d y  ¥  (HVC)

^I 

= 立 d x  ¥  ( ' HFx)    . .

^I 

+手 d y   (HFy)+S

山 )

ここで ，Hは全水深 ，Fx' Fyはx，y方 向 の 拡 散fluxを表す.また ，Dxx' D_XY' Dyx， Dyyを拡 散係数テンソルの各成分とすれば，Fx' Fyは一般的に次式で表せる.

a c  

^r. 

a c  

Fx  =  Dxx

^{て 一}

。

^X

+Dxy

て一^oy 

a c  

^r. 

a c  

Fy  =  Dyx

^{て 一}

+Dyy

τ‑

OX  oy 

(4.12a) 

(4.12b) 

4.3節におけるモデルは， (4.12)式では内湾における拡散能に等方性を仮定し ，Dxx 

= 

Dy_rry = D_‑， ， Dvv= Dvv _{‑Xy  ‑ } 宮}

= 

0とおき，場所毎の分散係数を (4.2)式で表した場合に相当する.

本節では，ここまで等方性が仮定されてきた分散係数に異方性を導入し，年平均的な 分散係数の異方性についての検証を試みる.

4.4.1非定常性の導入について

分散係数の非定常性を考慮するために，代表流速スケールとして各時刻の底面摩擦速

度U.(=九([p + v2 ^)ln ) を，代表長さスケールとして各時刻の全水深Hを導入 し，分 散係数を次式の様に表す.

D = 

r 

U*H  ^(4.13) 

博多湾における塩素イオン濃度の拡散計算により比例定数

r

のチューニングを行った ところ，最適値として

r

^= 3200が 得 ら れ た . 比 較 に 用 い た 実 測 デ ー タ は 福 岡 市 が 1987 年に湾内の 22ケ 所 に お い て 毎 月 定 点 測 定 し た デ ー タ を 基 に 年 平 均 値 を 算 出 し ， そ の 分 布に最もフィットする

r

を最適値と判断した.但し，計算手法は 4.3節と同じである が，移流項の計算スキームのみ風上差分を用いて有限差分法により計算した .

r  = 

3200  の場合の塩素イオン濃度分布の計算結果を図‑4.15a)に，実測値と計算値を比較したも のを図‑4.15b)に示す.

4.4.2異方性の導入について

移 流 分 散 効 果 は 流 れ 方 向 の み に 生 じ る は ず で あ る か ら ， 流 下 方 向

( c

‑方向) ，には 移流分散と乱流拡散による輸送を，流れに垂直な方向 (η‑方 向 ) に は 乱 流 拡 散 の み に よる輸送を考慮して拡散係数を次式の様に表した(図‑4.16参照) . 

vL  

ε 

+n ua  

a=  戸 め

D 

， ， 

(4.14a)  (4.14b) 

ここで ，DL， D Tは

c

，η方向の全拡散係数を表し， εlが移流分散係数を， εLt'εTtが

c

，η方 向の乱流拡散係数を表している.εl'ELt'εTt についても (4.13)式と同様なモデル化を行い 次式の様に表した.

f H H  

ι

﹂ 牢 牢

ι J U  

l L T  

k k k   一 一 一 一

一 一

向 ゐ 匂

(4.15a)  (4.l5b)  (4.15c) 

図‑4.16に示すように x軸と

C

軸がなす角を 0とすると， (4.12a， b)式中の Dxx等は 以下のように表せる.

Dxx  =  DL 

COS2

θ+DTsin

²

e 

Dy y 

= 

DL  s i n  

2e 

+  DT 

COS2e 

Dx  y  =  Dy  x  =  (  _DL  ‑DT )  s i n   e  ^{c o s} 。

(4.16a)  (4.16b)  (4.16c) 

このようなモデル化は以前にも Hollyand Usseglio‑Polatera (1982)やFalconerand Li u 

(1988)， Lin and Falconer (1997)などにより内湾の拡散シミュレーシヨンに対して行われ ていたが，前2者 は (4. 15 a， b， c)式 中 の 係 数k[，k_Lt， k_Ttとして Elder(1959)による開水路等 流場における値 k[= 5.86， kμ= 0.07， k_Tt = 0.23  (Falconer and Liuはk_Tt= 0.15 )を，後者は k，:::13， kμ= 0， k_Tt

= 

1.2を採用していた.

.博多湾に Hollyand Usseglio‑Polateraの 値 (k[ = 5.86， k_Lt = 0.07， k_Tt = 0.23 )を適用した場 合の計算結果を図‑4.17に示す.図‑4.15a)と比較しでも分かるように拡散能が弱すぎ るために湾奥への塩分の侵入が抑えられており，これらの係数の値が妥当でないことが 分かる.これは，移流分散係数の与え方が小さ過ぎたためと考えられるので，係数 k[の 値のみを変化させることにより博多湾の実測塩素イオン濃度分布の再現を試みた.最適 値として k[=2000が得られたが，その場合の計算結果(図‑4.18)と現況を良く再現し ている図‑4.15の比較から明らかなように，湾奥への塩分の侵入は起こっているもの の，非等方性が強すぎたため濃度分布はかなり歪んだものとなっている.また，k[ を 2000より更に大きくした場合についても計算を試みたが，これ以上非等方性が強くなる

と数値振動が発生してしまい不合理な計算結果となってしまった.これについては未だ に原因が不明であり，検討の余地が残されている.

4.4.3異方性の修正について

4.4.2のモデルでは非等方性が強すぎたため非現実的な濃度分布が計算された.そこ で，分散係数の非等方性の度合いを変化させることにより修正を試みた.修正のための 係数

α $

を導入し，D_L， DTを次式の様に表す.

D

_L =

ε ' [   + ε

Lt  D T =

α

^牢

ε ' [+ ε

Tt 

e ~

(4.17a)  (4.17b) 

α. 

= 

0が 完 全 な 非 等 方 性 ， 氏 =1が等方的な場合に相当する.そこで，k[を等方的な 場合の最適値 k[= 3200に固定し，

α .

を 0から 1まで変化させて濃度分布の計算結果の 比較を行った.比較を容易にするために，福岡市の測定点22ケ所で計算結果と実測結果 の聞の絶対誤差の平均値をとった(表

‑4.3

参 照 ) と こ ろ ， 等 方 的 で あ る 久 =

1

の場合 が最も精度が良いという計算結果が得られた.これより，年平均の分散係数はかなり等 方的であることが確認された.

ドキュメント内 平面的流れ場・拡散場の高精度計算法の開発 (ページ 65-76)