森林のＣＯ２フラックス測定法
【課題】 森林全体におけるＣＯ２の吸収量を高精度に評価する。
【解決手段】 シンチロメータを用いた場合と渦相関法による場合の温度変動の消散率の比εＴSAS／εＴECが、シンチロメータを用いた場合の測定範囲ＡSASと渦相関センサによる場合の測定範囲ＡECの比ＲＡ（＝ＡSAS／ＡEC）との相関、および大気の乱流強度σｗ／ｕとの相関が高いことに基づき、これら２つの要素ＲＡおよびσｗ／ｕをパラメータとする関係式をあらかじめ用意しておくとともに、大気の乱流強度σｗ／ｕを測定し、渦相関法によって求めた顕熱フラックスＨECをシンチロメータによって求めた場合に相当する顕熱フラックスＨSASに補正し、尚かつ、渦相関法によって求めた潜熱フラックスλＥECおよびＣＯ２フラックスＦｃECに関しても、シンチロメータによって求めた場合に相当する潜熱フラックスλＥSASおよびＣＯ２フラックスＦｃSASに補正するようにした。
【解決手段】 シンチロメータを用いた場合と渦相関法による場合の温度変動の消散率の比εＴSAS／εＴECが、シンチロメータを用いた場合の測定範囲ＡSASと渦相関センサによる場合の測定範囲ＡECの比ＲＡ（＝ＡSAS／ＡEC）との相関、および大気の乱流強度σｗ／ｕとの相関が高いことに基づき、これら２つの要素ＲＡおよびσｗ／ｕをパラメータとする関係式をあらかじめ用意しておくとともに、大気の乱流強度σｗ／ｕを測定し、渦相関法によって求めた顕熱フラックスＨECをシンチロメータによって求めた場合に相当する顕熱フラックスＨSASに補正し、尚かつ、渦相関法によって求めた潜熱フラックスλＥECおよびＣＯ２フラックスＦｃECに関しても、シンチロメータによって求めた場合に相当する潜熱フラックスλＥSASおよびＣＯ２フラックスＦｃSASに補正するようにした。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、森林のＣＯ２フラックス測定法に関する。さらに詳述すると、本発明は、渦相関法（ＥＣ）により評価されていた森林全体のＣＯ２吸収量の評価技術の改良に関する。
【背景技術】
【０００２】
温暖化ガス削減策として京都メカニズムにおける森林吸収源を利用するには、森林吸収源におけるＣＯ２の吸収量の評価方法を示し、実行する必要がある。温暖化抑制や気候変動における森林の影響を定量化する上では、森林全体のＣＯ２収支の評価（フルカーボンアカウンティング）が重要である。このような森林全体のＣＯ２収支は、ＣＯ２の交換量（フラックス）を積算することで評価できるものであり、現在、地球規模のＣＯ２収支の評価・予測が盛んに行われている。このように評価や予測を行うにあたっては、地表面で実測されたエネルギー・物質交換量（フラックス）の検証データとしての利用価値が高いことから、フラックスサイトをネットワーク化し、広域的かつ長期間にわたってフラックスデータをグラウンドトゥルースとして利用できる環境の整備が世界的に進められている。また、現在、ＩＰＣＣ（Intergovernmental Panel on Climatic Change、1988年に設立された国連の組織で、地球温暖化問題についての科学的な議論を行う場）が推奨するＣＯ２吸収量の評価手法があるが、この評価手法は森林のＣＯ２収支のうち土壌からのＣＯ２放出を正確に把握できないことから、森林のＣＯ２収支を正しく評価できる新しい手法が今後は求められるであろうという面もある。
【０００３】
例えば、小規模な公園程度の大きさであればＣＯ２の出入りの量に差が生じるが、数十ｈａあるいは数百ｈａといった大規模な森林であれば出入りの量はほぼ等しいと考えることができる。したがって、森林におけるＣＯ２吸収量を評価するにあたり、森林が横方向に均一だと考えれば、森林に対して横方向から流入する風と横方向に流出していく風のＣＯ２の流量（交換量）はほぼ等しいとして扱うことが可能となる（図１参照）。したがって、この考えに基づけば、横方向の流量は考慮せず、森林の上部空間において上下方向に出入りする空気量などを例えば１年間測定することにより、１年の間に森林がどの程度ＣＯ２を吸収したか、あるいは放出したかがわかる（図１参照）。より具体的には、森林の上下方向におけるＣＯ２の交換量速度を把握し続け、最後に積算すればある期間（例えば１年）の吸収量（放出量）がわかる。この吸収ないしは放出の量、つまりはＣＯ２の交換量を単位時間あたりの量として考えたものが「フラックス」、つまりいわば交換速度ということになる（図１参照）。
【０００４】
ここで、森林上におけるＣＯ２の収支を表すＮＥＥ（Net Ecosystem Exchange：生態系純交換量）は、ＣＯ２フラックスＦＣを用いて次式で表される。
［数７］
ＮＥＥ＝ＦＣ＋ＦＳ
ＦＳは森林内の空間に一時的に存在するＣＯ２量（林冠内貯留量）で、１日の中で貯留と放出を繰り返し、１日を通した収支は非常に小さい。したがって、１年間のＮＥＥは、ＣＯ２フラックスの時系列１年分を積算することで求めることができる。１年間のＮＥＥ（mgＣＯ２m-2yr-1）は光合成の同化作用によって森林生態系が１年間に吸収したＣＯ２量を示す（図５参照）。
【０００５】
従来、このＣＯ２フラックスを測定するための手法の一つとして、渦相関法（ＥＣ：Eddy-covariance method）が広く用いられている。この渦相関法は、垂直風速の変動（ｗ’）と気温や湿度、ＣＯ２密度変動の測定値（ｃ’）の共分散の平均値（ｗ’ｃ’（文中ではアッパーラインを付していないが、算術平均を表しているものとする））から乱流輸送量であるフラックスを求める、つまり、図１に示してあるように、
ＣＯ２フラックス＝風の乱れ×ＣＯ２密度変動＝ｗ’ｃ’
という式に従い、ＣＯ２フラックスを求めるというものである。この渦相関法による測定手法は、各物理量を直接測定するものであることから最も仮定が少ない測定法とされている。また、この渦相関法を利用する場合の施設としては森林上空の風の流れとＣＯ２の密度の変動を直接的に測定するための施設、例えば森林内に建てられたタワーや当該タワーに設けられたセンサなどが利用されている（例えば、非特許文献１参照）。センサとしては、風の乱れを測定するための超音波風速計や、ＣＯ２密度変動を測定するためのＣＯ２／Ｈ２Ｏ分析計などが用いられる。
【０００６】
【非特許文献１】Thiermann, V. "A displaced-beam scintillometer for line-averaged measurements of surface layer turbulence". Tenth symposium of turbulence and diffusion, 29 Sept -2 Oct, 1992, Portland, OR., published by the American Meteorological Society, Boston, MA.: ｐ244-ｐ247 (1992)
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、従来の測定法の場合、測定精度が十分ではないとの見方がある。すなわち、フラックスは大気下層の地表面付近では大気乱流によって駆動されることから、晴天日の夜間など大気が安定した状態や、大気の対流が大きい不安定な条件では、渦相関法による熱フラックスと放射計による放射フラックスが釣り合わない現象（インバランス）が多発している。そして、このような現象（インバランス）は、主に渦相関法に起因する問題と認識されている。
【０００８】
そこで、本発明は、森林全体におけるＣＯ２の吸収量を高精度に評価するための森林のＣＯ２フラックス測定法を提供することを目的とする。
【課題を解決するための手段】
【０００９】
かかる目的を達成するため、本発明者は種々の検討を行った。森林全体におけるＣＯ２の吸収量を高精度に評価するため、森林全体における時間平均ないしは空間平均のＣＯ２フラックスを測定することを考えた。ただし、上述したように、渦相関法だと測定精度が十分ではないとの見方があり、本発明者はこの問題について検討した結果、森林上空の１点にて測定を行ういわば「点測定」型であることに理由があるのではないかと考え、さらに検討を進めた結果、パス（トランスミッタとレシーバとの設置間隔）を長くとれるシンチロメータを使ってＣＯ２フラックスをいわば面的に広い領域にて測定するという着想を得るに至った。シンチロメータならば、トランスミッタとレシーバとを数十ｍあるいは100ｍ以上離して設置し、それらの間におけるレーザ光の屈折率変化を測定することにより、少なくとも渦相関法の場合よりも広いエリアについてＣＯ２フラックスを測定することが可能となる。
【００１０】
ただし、シンチロメータの場合だと、温度の乱流強度しか測定することができないために、上述したＣＯ２フラックスを求めるための要素である「ＣＯ２密度変動」に関するデータを直接求めることができないという欠点もある。つまり、より広い領域にてフラックス測定することにより測定精度を高めることが期待できても、ＣＯ２フラックスを求めるための要素を直接求めることができなければそもそもＣＯ２フラックス測定が不可能となってしまう。
【００１１】
ここで、本発明者は、パスを長くとれるシンチロメータを測定手段として利用しつつＣＯ２フラックスを求めるためさらなる検討を加えた。そこで、測定手法の差異に対するソースエリアの影響を調べるために、ソースエリアおよびシンチロメータと渦相関センサによる消散率を比較することに想到し、さらなる検討の結果、ある知見を得るに至った。
【００１２】
本発明はかかる知見に基づくものであり、まず、森林のＣＯ２フラックス測定法として、森林上空のシンチロメータのトランスミッタから照射したレーザ光をシンチロメータのレシーバで受光する際、森林上における乱流あるいは温度変動に起因するシンチレーションに基づく光強度の変動を計測して温度変動の消散率εＴを求め、さらに、森林上空の渦相関センサによって、温度に関する変動の分散の平方根σＴ、水蒸気に関する変動の分散の平方根σｑ、ＣＯ２に関する変動の分散の平方根σｃを測定し、その後、以下の関係式
【数８】
【数９】
に基づき、水蒸気の消散率εｑ、ＣＯ２の消散率εｃをそれぞれ算出し、スカラー量である各消散率εＴ、εｑ、εｃからフラックスを求めることを特徴とするものを発明した。
【００１３】
例えば森林の上空で空気乱流や熱変動が生じると、光の散乱が生じ、光強度が変動する。そこで、このように相関的に変動した光強度に基づき、乱流運動エネルギーの消散率、あるいは温度変動の消散率を求める（図２参照）。また、空気の乱れとＣＯ２の密度とを掛け合わせる、あるいは空気の乱れと温度変動と掛け合わせることによって熱のフラックスを求めることができる。ここで、フラックスの生成と消散には一対の関係がある。そこで、消散率を求めることによって間接的にフラックスを算出することが可能である。
【００１４】
ただし、この手法においてはＣＯ２のフラックス自体を直接的に求めることはできず、顕熱フラックス（伝熱）と風速シアー（乱流の強さ）という２つが直接的に求められるに過ぎない。そこで、この方法を応用してＣＯ２のフラックスあるいは水蒸気のフラックスを求める手法を新たに導入する。すなわち、上述したようにフラックスの生成と消散とには一対の関係があることから、本発明においては、消散率を求めることによってフラックスを間接的に算出することとしている。この場合、渦相関センサによる計測値からシンチロメータの空間スケールに対応するフラックスを算出する点が特徴的である。
【００１５】
より具体的に説明すると、シンチロメータを用いればフラックスを間接的に測定することは可能だが、このシンチロメータの場合、温度と乱流の強さしかわからない。そこで、温度と水蒸気の消散率の比が、温度と水蒸気の分散の比の二乗と等しいという関係式を導入し、渦相関センサによって計測したデータから、温度に関する変動の分散の平方根σＴと水蒸気に関する変動の分散の平方根σｑの比を求め、さらには温度に関する変動の分散の平方根σＴとＣＯ２に関する変動の分散の平方根σｃの比を求め、ここから、シンチロメータによっては直接的には測りえない水蒸気の消散率εｑ、ＣＯ２の消散率εｃをそれぞれ算出することとしている（図３参照）。こうしてスカラー量である各消散率εＴ、εｑ、εｃを求めることができれば、これらの結果からフラックス（例えばＣＯ２フラックス）を求めることが可能となる。
【００１６】
なお、渦相関法において熱フラックスの過小評価が見られたことがあった。このような測定手法に対し、本発明の場合、シンチロメータによる熱フラックスは熱収支のインバランスを緩和させ、より妥当な評価結果が得られる。また、ＣＯ２フラックスは顕熱および潜熱フラックスと同じ測定原理に基づくことから、シンチロメータの適用によって、より信頼できる値が得られることになる。
【００１７】
ここで、本発明者は検討を重ねた結果、森林上に適用したシンチロメータから得られた乱流運動エネルギーおよびスカラー分散の消散率（フラックス算出のための支配的パラメータ）が、測定パスに直交する風向時に渦相関センサによる消散率を特徴的に上回る現象を知見するとともに、この現象の原因のひとつに、シンチロメータと渦相関センサの間でソースエリアが大きく異なることがあると考えた。かかる知見に基づいて、上記のようなＣＯ２フラックス測定を実施する場合には、シンチロメータにおけるレーザ光の照射方向たるパスの方向が森林における風向と直交するようにトランスミッタおよびレシーバを配置しまたは配置を変換することが好ましい。シンチロメータによる測定領域を広げるという観点からは、上記のようにパス方向を風向に直交させれば、シンチロメータを利用してＣＯ２フラックスを測定する場合により大きな効果を得ることが可能となる（図４参照）。例えば、風向が一定ないしは安定となっている場合にパスを直交させることとすれば測定精度が高くなる。また、特に測定パスをより長距離化した場合に、シンチロメータは広い領域においてフラックスが導出できるため、より広い面的な測定が可能となるという利点がある。本発明の場合、トランスミッタとレシーバの少なくとも一方を移動可能としておけば風向に応じてパスの方向を適宜変化させることが可能となる。
【００１８】
以上がシンチロメータを利用した森林のＣＯ２フラックス測定法であるが、本発明者はさらに検討を重ねた。シンチロメータを利用すれば上記のような効果が得られる反面、レーザ光の光路が遮られると測定不能となるため、雨天や霧など悪天候時のデータが得られないという面もある。また、シンチロメータの測定パスに風向が直交する際には領域を代表する値が得られる一方で、パスに平行な風向時には渦相関法と同等の評価値となる。フラックスの積算によってＣＯ２吸収量を得るためには、シンチロメータの空間スケールを持つフラックスを長期に安定して算出する必要もあると考えた。そこでさらに検討を重ねた結果、シンチロメータと渦相関法による温度変動の消散率の比εＴSAS／εＴECが、(i)センサの測定範囲の比（換言すればソースエリアの比、ＲＡ＝ＡSAS／ＡEC ）および (ii)大気の乱流強度（σｗ／ｕ） と相関が高いことから、これら２つの要素をパラメータとして、渦相関法の測定値をシンチロメータ相当のフラックスに変換する手法に想到するに至った。
【００１９】
請求項１に記載の発明はこのような考えに基づき想到したものであり、森林全体のＣＯ２吸収量を測定するためのＣＯ２フラックス測定法において、シンチロメータを用いた場合と渦相関法による場合の温度変動の消散率の比εＴSAS／εＴECが、シンチロメータを用いた場合の測定範囲ＡSASと渦相関センサによる場合の測定範囲ＡEC の比ＲＡ（＝ＡSAS／ＡEC ）との相関、および大気の乱流強度σｗ／ｕ との相関が高いことに基づき、これら２つの要素ＲＡ（＝ＡSAS／ＡEC ）およびσｗ／ｕ をパラメータとする関係式
【数１０】
を、シンチロメータおよび渦相関センサを使って実際に得られた一定期間における測定結果に基づきあらかじめ用意しておくとともに、数式１０中における測定範囲の比ＲＡには一定値を与えることによりすべての風向に対するシンチロメータの広い測定範囲を適用することとし、その後、森林の上空に実際に設置した渦相関センサによって大気の乱流強度σｗ／ｕを測定し、以下の関係式
【数１１】
に基づき、渦相関法によって求めた顕熱フラックスＨECをシンチロメータによって求めた場合に相当する顕熱フラックスＨSASに補正し、尚かつ、渦相関法によって求めた潜熱フラックスλＥECおよびＣＯ２フラックスＦｃECに関しても、以下の関係式
【数１２】
【数１３】
に基づき、シンチロメータによって求めた場合に相当する潜熱フラックスλＥSASおよびＣＯ２フラックスＦｃSASに補正し、シンチロメータにて測定したのと同様の測定結果を得ることを特徴としている。この場合のＲＡは風向によって変化するシンチロメータの空間スケール、σｗ／ｕは大気の乱流状態を表している。
【００２０】
本発明の場合、上述の(i)と(ii)を変数とする関数（数式１０）を、実際にシンチロメータおよび渦相関センサを使った過去の測定結果に基づきいわば経験値としてあらかじめ用意しておく。その後の計測は渦相関法のみを実施し、シンチロメータを使った実測はもはや行わないが、あらかじめ用意された関係式を利用し、フラックスを補正することによってシンチロメータを使った場合と同様の測定結果を得ることができる。つまりは、当該測定対象森林にシンチロメータを実際に設置する必要はなく、あくまでも従来と同様の手順により渦相関センサのみにより実際の測定を実施して測定データを得ることができる。例えば渦相関センサのみによる測定結果は実際の値よりも過小評価される傾向があるが、本発明においては特定の関係式に基づきこの結果を補正ないしは補間し、あたかもシンチロメータを使って得たのと同様、真の値により近い測定結果を算出することができる。
【発明の効果】
【００２１】
前述のＣＯ２フラックス測定法によれば、測定パスを長距離化できるシンチロメータを利用することにより、広い領域においてフラックスを計測することが可能となる。例えば従来の渦相関法の場合、いわば点で測定するため測定領域が極めて狭いのに対し、本発明にかかるフラックス測定法の場合だと大きい領域を測定できるため、狭小領域における変動には依存しない平均的かつ精度の高い測定データを得ることができる。これにより、森林全体のＣＯ２吸収量を高精度に評価することが可能となる。
【００２２】
また、前述のＣＯ２フラックス測定法においては、風向にパスを直交させるか、あるいはこれに近似した状態とすることによってより広い測定領域を確保することができる（図４参照）。これによれば、森林全体のＣＯ２吸収量をより高精度に評価することが可能となる。
【００２３】
請求項１の発明によると、シンチロメータは使わずに、シンチロメータを使ったのと同様の高精度の測定結果を得ることができる。こうした場合、トランスミッタとレシーバのそれぞれを設置するためのタワーを設けなくても済むため、コストを抑えつつ、ＣＯ２吸収量を高精度に評価することが可能となる。要は、測定領域の異なる計測値から導いた補正を施すことによって、あたかもシンチロメータにて計測を実施したのと同様の結果を得ることができ、渦相関法による熱フラックスのインバランスを解消に向かわせることが可能となる。
【００２４】
しかも、シンチロメータは光学的な計測機器であるために雨や霧の場合に稼働率が下がるといったように天候に左右されやすいが、実際にシンチロメータを使わずに高精度な結果が得られるようにした本発明にかかる測定法によれば天候の影響を受けずに随時測定することが可能になるという利点もある。
【発明を実施するための最良の形態】
【００２５】
以下、本発明の構成を図面に示す実施の形態に基づいて詳細に説明する。
【００２６】
図１〜図８に本発明の実施の一形態を示す。本発明にかかるＣＯ２測定法は、森林全体におけるＣＯ２フラックスを計測し評価するというもので、森林上空のシンチロメータ１のトランスミッタ２から照射したレーザ光をシンチロメータ１のレシーバ３で受光する際、森林上における乱流あるいは温度変動に起因するシンチレーションに基づく光強度の変動を計測して温度変動の消散率εＴを求めるか、あるいは経験値としてあらかじめ用意しておき、さらに、森林上空の渦相関センサによって、温度に関する変動の分散（またはその平方根σＴ）、水蒸気に関する変動の分散（またはその平方根σｑ）、ＣＯ２に関する変動の分散（またはその平方根σｃ）を測定し、その後、以下の関係式
【数１４】
【数１５】
に基づき、水蒸気の消散率εｑ、ＣＯ２の消散率εｃをそれぞれ算出し、スカラー量である各消散率εＴ、εｑ、εｃからフラックスを求めるというものである。
【００２７】
まず、ここで本明細書中における各種用語について説明しておく。まず、一般に「フラックス」といった場合、気流の乱れによって鉛直方向に運ばれる乱れのエネルギー（運動量フラックス）や、熱・物質の単位面積あたりの移動速度（顕熱フラックス、潜熱フラックス、ＣＯ２フラックス）のことを意味する。本明細書では、単位時間当たりのＣＯ２などの交換量のことをフラックスと呼ぶ。「渦相関法」は、超音波風速計とＣＯ２／Ｈ２Ｏ変動計を用いて計測した水平および鉛直方向の風速と物質量の変動成分からフラックスを算定する手法である。「消散率」は、大気の乱れのエネルギーや、熱・物質の密度変動の消散を表す値である。接地境界層の乱流エネルギーについての生成と消散に関する収支式から乱流フラックスが算定される。
【００２８】
本実施形態では、まずシンチロメータ１によるフラックスの算出方法の概要について説明する。次に、その手法を拡張し、シンチロメータ１の空間スケールから決定した関数を用いて渦相関法の測定値を補正するというフラックス算出方法について説明する。
【００２９】
シンチロメータ１（本実施形態においてはSAS：Small Aperture Scintillometerを用いている）は、可視レーザーを射出するトランスミッタ２とレーザー強度を検出するレシーバ３の一対で構成されており、乱流によって引き起こされる大気の屈折率の変動を光強度変化として捉え、顕熱フラックス算出のパラメータである温度変動の消散率（εＴ）を算出する（図２参照）。温度変動の消散率は、渦相関センサによって計測した風速と温度の時系列からも計算することができる。フラックスの収支式にモニンオブコフの相似則を適用することで、消散率は次式によってスカラーフラックスである顕熱Ｈ（Ｗｍ−２）、潜熱λＥ（Ｗｍ−２）およびＣＯ２フラックスＦＣ（mgＣＯ２m−２s−１）と関連付けられる。
【数１６】
【数１７】
【数１８】
ここで、ｕ＊は摩擦速度（ｍｓ−１）、ρは空気密度（ｇｍ−３）、ｋｋａｒはカルマン定数（＝0.4）、ｚは測定高さ（ｍ）、ｃｐは空気の定圧比熱（Ｊｇ−１Ｋ−１）およびλは水の蒸発潜熱（Ｊｇ−１）を表す。εｑ、εｃはそれぞれ、比湿（ｑ ｋｇｋｇ−１）変動およびＣＯ２密度（ｍｇＣＯ２ｍ−３）変動の消散率である。φε(ζ)は無次元化された消散率で大気安定度を表すパラメータ
［数１９］
ζ＝ｚ／Ｌｍｏ（Ｌｍｏ：モニンオブコフ長）
の関数として与えられる。
【００３０】
続いて、シンチロメータ１によって決定したｕ＊およびεＴと、同時に測定したｑおよびＣＯ２密度変動の分散から、εｑ，εｃを求めて潜熱フラックスおよびＣＯ２フラックスを算出する手法の内容について説明する。
【００３１】
シンチロメータ１によるフラックス算定は消散法に基づいて行われる。乱れは生成消散によって成り立っており、渦相関法がスカラー変動や乱流運動エネルギー（気流の乱れのエネルギー：Turbulent Kinetic Energy＝ＴＫＥ）生成の直接測定に相当するのに対し、消散法はスカラー変動やＴＫＥの消散率（エネルギーの消散を表す正の値）からそれらの生成を計算する手法である。消散法は船の動揺など低周波ノイズの影響を受けにくいことから海上でのフラックス測定などに応用されてきた。消散法によるフラックス算出では、まずＴＫＥとスカラーの変動の消散率を求め、次いでフラックスと消散率との収支式からフラックスを算定する。シンチロメータ１はこのうちＴＫＥとスカラー変動の消散率を測定する装置である。
【００３２】
消散法の原理について以下に説明する。水平一様で定常な接地境界層における乱流運動エネルギーＴＫＥ(ｅｔ)およびスカラー変動の収支方程式は以下のように書かれる。
【数２０】
【数２１】
【００３３】
ここでｓは気温Ｔ、比湿ｑ、やＣＯ２などのスカラーを表し、ｐは気圧で
［数２２］
θ’ｖ＝θ（1＋0.61ｑ）
は仮温位（ただしθは温位）、ｅｔはＴＫＥ（ｅｔ＝（ｕ２＋ｖ２＋ｗ２）／２）を表す。ここで、ｕ，ｖ，ｗは平均風向、横風方向、鉛直方向の風速成分である。また、この文中ではアッパーラインを付していないが、数式中のバーは平均値、プライムは平均値からの変動を表す。数式２０の各項はそれぞれ生成項、浮力項、輸送項、圧力項、そして消散率である。また、数式２１の各項は、生成項、輸送項、そして消散率である。両式第１項のｕ’ｗ’およびｗ’ｓ’（ただし本文中ではアッパーラインを省略）が求めるべき運動量およびスカラーのフラックスに相当する。接地境界層においては両式とも輸送項および圧力項を無視できると仮定して、数式２３で定義される摩擦速度ｕ＊、
【数２３】
Monin-Obukov長（Lmo）およびKarman定数（ｋkar 0.4）を用いて、それぞれを
ｕ＊3／ｋkar(ｚ−ｄ0) と ｋkar(ｚ−ｄ0)ｕ＊／（ｗ’ｓ’）２ によって無次元化すると、数式２０、数式２１の各項は測定高さｚ(m)において、
【数２４】
【数２５】
と表される。φεとφεｓはＴＫＥおよびスカラーの無次元消散率を、φｍとφｓはＴＫＥおよびスカラーの無次元した鉛直勾配を表す。ｄ0 はゼロ面変位で、植生地上での風速プロファイルの対数分布を仮定した際の、高さを測る基準面を上方に修正するためのパラメータである。群落高×0.7という値がよく使われる。さらに、ここでＴＫＥおよびスカラー変動の消散率は次のように定義される（Brutsaert，1982）。
【００３４】
【数２６】
【数２７】
水平一様な接地境界層ではこれらの無次元消散率は、大気安定度を示すパラメータζ＝（ｚ−ｄ０）／Ｌｍｏ（安定時ζ＞０−、不安定時ζ＜０）の関数として表されることを数式２４、数式２５は示している。ここで用いられるMonin-Obukov長Ｌｍｏは、以下の式で定義される。
【００３５】
【数２８】
【数２９】
ρは空気密度、Ｔｖは仮温度であり、Ｈ、λＥは地表面における顕熱および潜熱フラックスである。数式２６と数式２７のｓ＝θおよびｓ＝ｑに関する式をそれぞれｕ＊、Ｈ、λＥについて解くと、
【００３６】
【数３０】
【数３１】
【数３２】
で、数式３０〜数式３２は、無次元消散率の関数形が既知であれば、消散率を与えることで連立方程式として解くことができ、フラックスｕ＊、Ｈ、λＥが求められる。大気安定度パラメータζによる無次元関数表現は、接地境界層におけるMonin-Obukov similarity（ＭＯＳ）が成立することに基づいており、物理量の鉛直プロファイルがζで一意に表現できることを示している。例えば数式２４、数式２５中のφｍ≡ｋkarｚ／ｕ＊∂ｕ／∂ｚ、
φｓ≡ｋkarｚｕ＊／ｗ’ｓ’∂ｓ／∂ｚ をζの関数で表すのにBusinger-Deyer式がよく用いられる。また、無次元消散率の関数形も観測結果を元にいくつか提案されている。例えば、Kader(1992)による経験式は以下である。
【００３７】
【数３３】
【数３４】
温位と比湿の大気安定度に対する応答が相似であれば φεｃ＝φε０＝φεｑ が成立し、数式３０、数式３１、数式３２は繰り返し計算によって解くことが出来る。更にスカラーの無次元消散率の相似性の成立を仮定した場合は、ＣＯ２の消散率から、数式３１、数式３２と同様にしてフラックスを算出することが出来る。
【００３８】
次に、渦相関センサによる消散率の算定について説明する。シンチロメータ１との比較に用いるために、渦相関センサによる計測値から消散率を算定する方法について述べる。消散率の算定には、スペクトル（エネルギーや物質量変動の周波数表現）の慣性小領域における乱流特性を利用して間接的に求める慣性消散法（inertial dissipation method）を用いる。ただし、フーリエ変換の結果は一般に乱れが大きく、消散率の推定誤差要因となりうるため、ここではスペクトルの実空間表現である構造関数にもとづいた消散率の推定法を用いる。ＴＫＥおよびスカラー変動のスペクトルに対応する２次の構造関数Ｄab(r)は次式で定義される。
【数３５】
ここでａ，ｂはそれぞれｕ，θ，ｑなどの風速やスカラーを表し、ｒは実空間における距離である。時系列のデータ間の時間差をτとすると平均風速を用いて、τｕ＝ｒの関係にある。Kolmogorovの理論によればｒが慣性小領域にあるとき２次の構造関数と消散率は以下の関係にある。
【数３６】
【数３７】
ここでεとεｓはそれぞれＴＫＥおよびスカラー変動の消散率を表す。αuuはKolmogorov定数（およそ0.55）で、αssはObukhov-Corrsin定数（およそ0.8）と呼ばれる普遍的な定数である。比例定数Ｃs２は構造関数定数で、次のａとｂの組み合わせ［Ｔ：Ｔ］（気温）、［ｑ：ｑ］（比湿）および［Ｔ：ｑ］はそれぞれ構造関数定数ＣＴ２、Ｃｑ２およびＣＴｑに対応する。数式３６、数式３７に示すように２次の構造関数がｒ2/3に比例していることから、慣性小領域に相当するｒの範囲で構造関数にｒ2/3の曲線をフィットさせることで数式３６、数式３７から消散率を求め、消散法によるフラックス計算へと進む（図１１参照）。超音波風速計のパス長ｄｓより小さい乱流変動は測定されず、測定高さよりも大きい渦はその一部分しか測られないことから、十分に乱流変動計測が保障されるフィット範囲として ２ｄｓ＜ｒ＜（ｚ−ｄ０）／２ を選択する。
【００３９】
次に、シンチロメータ１を用いた消散率の算定方法について述べる。SLS40A（Scintec AG，Germany）は、市販のdisplaced-beam small aperture scintillometer で、波長0.67μｍのHeNeレーザーを偏光面が直交する２本のビームに分離し、２光路のシンチレーションからＴＫＥおよび気温変動の消散率を求め、フラックスを算出する装置である。SLSシリーズは単一光源を用いることで動作および計測の安定性と経済性を両立させた製品である。乱流による屈折率変動は大気中を伝搬する単色光の光強度のlog分散として計測することが出来る。独立した２光源と受光機がある場合の光強度のlog共分散Ｂ12は次式で表される（Thiemann，1992）。
【００４０】
【数３８】
ここで、ｘは、パス長Ｌの光軸上での位置；Ｋ＝２π／λは光の波数；ｋ乱流の波数；Ｃｎ２は屈折率変動の構造関数定数、ｌ０内部スケール（乱流渦が徐々に小さな渦に崩壊していく際、これ以上小さくなると分子粘性により渦が維持できなくなり熱として分解していくスケールの目安）である。Ｊ０およびＪ１は０次と１次の第１種ベッセル関数で、ビーム間の距離ｄおよび受光機の直径Ｄによって定まる感度分布を表す。検出器単体による分散を表すＢ１とＢ２は数式３８でｄ＝０とすることで得られる。数式３８中のφｎは屈折率変動の次元スペクトルで、
【数３９】
であり、ｆ（ｋｌ０）は、Hill and Clifford(1978)による屈折率変動の３次元スペクトルモデルである。数式３８は、σｘ２＜0.3の低散乱条件（Rytov approximation）が前提である。式中の被積分関数はシンチロメータ１の感度分布を表す荷重関数で、パスの中央で最大値となる釣鐘型の関数形を取る。数式３９を数式３８に代入し、ｘおよびｋについて積分することで次式が得られる。
【００４１】
【数４０】
ｆＢはｌ０，ｄおよびＤの増加に対するＢ12の減衰を表す関数である。相関係数ｒ12＝Ｂ12／（Ｂ１Ｂ２）１／２ ＝Ｂ12／Ｂ1＝Ｂ12／Ｂ2 はｌ０，ｄおよびＤの関数となり、ｄとＤを定数と考えると、相関係数ｒ12から内部スケールｌ０が求められる。ｌ０がわかれば、Ｂ1（あるいはＢ２）からＣｎ２が求められる。Scintecのソフトウェアでは、ｌ０に対するｒおよびｆＢをパス長Ｌに対する数表として用意し、Ｂ12とＢ1およびＢ2 からリアルタイムでｌ０とＣｎ２を計算する。得られたｌ０から乱流運動エネルギー（ＴＫＥ）の消散率は次の関係によって求められる。
【００４２】
【数４１】
ここでν（ｍ2ｓ-1）は空気の動粘性係数である。屈折率の変動に最も寄与するのは気温Ｔ（Ｋ）であり、ついで比湿ｑ（kgkg-1）の変動である。ＣＯ２などのガスは大気中に微量しか存在しないため屈折率変動には影響しない。屈折率変動の構造関数定数Ｃｎ２は気温変動の構造関数定数ＣＴ２（Ｋ２ｍ−２／３）、比湿変動Ｃｑ２（ｍ−２／３）および気温と比湿Ｔ−ｑの結合構造関数定数ＣＴｑ（Ｋｍ−２／３）によって次のように表される（Hill et al.，1980）。
【００４３】
【数４２】
定数ＡＴとＡｑは光の波長、気圧、気温および比湿の影響を受ける。本実施形態で使用するシンチロメータ１で利用する波長λ＝0.67μｍにおいては、Andreas(1989)による以下の式を用いる。
【数４３】
ここでＰ（ｈＰａ）は気圧である。光源に単一波長を用いるシンチロメータ１ではＣｎ２の値が一つしか得られないため数式４２からＣＴ２およびＣｑ２を一意に求めることが出来ない。この問題に対し、Thiermannは比湿がＣｎ２に与える影響が気温に比べて１オーダー小さいためこれを無視できるとし、数式４２の右辺第１項のみからＣＴ２を求めた。本実施形態で用いるシンチロメータ１に付属のフラックス計算ソフトウェアはこの計算方法を採用しており、運動量フラックスおよび顕熱フラックスのみが計算される。
【００４４】
また、本実施形態では水蒸気やＣＯ２密度などのスカラーについてもシンチロメータ１に対応する値を得るための手法を適用する。
【００４５】
まず水蒸気補正のための既存手法であるが、数式４２から明らかなように、シンチロメータ１は水蒸気の変動に対しても感度を持ち、水蒸気フラックスが顕熱フラックスに比べて非常に大きい条件ではＣＴ２が過大評価されてしまうことから水蒸気に対する補正の必要性が指摘されている（Green and Hayashi，1998）。水蒸気影響の古典的な補正法はWesely(1976)によって提唱されたボーエン比（顕熱と潜熱の比：Ｈ／λＥ）補正である。ボーエン比補正にはｗ’Ｔ’とｗ’ｑ’が必要なため、渦相関システムなどのセンサによる測定値を用いる。Moene(2003)はボーエン比補正について精査し、気温と比湿の標準偏差σＴ，σｑの比を用いる手法が、フラックスの比を用いた手法に対して、測定項目が少ないため誤差が小さいことを示している。すなわち
【数４４】
とすることで、数式４２は、Ｃｎ２とＣＴ２に関して次のように変形できる。
【数４５】
Ｔとｑの平均値と分散、そして相関係数ＲＴｑを他のセンサで計測することで、水蒸気の影響をキャンセルしたＣＴ２を計算することが出来る。
【００４６】
次にスカラーフラックスへの拡張であるが、水蒸気補正の際に導入した水蒸気と気温の構造関数定数比と分散比の相似性が成立するなら数式４４にしたがって、シンチロメータ１によるＣＴ２と、渦相関センサによるσＴ，σｑからＣｑ２を算出できると考えられる。同様にこの相似性をＣＯ２などのガスフラックスまで拡張すると以下が成り立つ。
【数４６】
これにより、渦相関センサによって得られた気温、水蒸気密度およびＣＯ２密度の分散とシンチロメータ１によって得られた気温変動の構造関数定数から水蒸気およびＣＯ２変動の構造関数定数ＣＴｑ２，Ｃｃ２が求められる。この算出法は純粋なシンチロメータ１による導出ではなく、渦相関センサとの組み合わせによることを留意する必要がある。しかし、現時点でシンチロメータ１と同様のスケールでの変動を測定する手法は存在しない。数式３７にｑおよびｃに関する構造関数定数を代入して消散率εｑとεｃが求められ、ε，εＴを加えて消散法によるフラックス計算へと進む。
【００４７】
続いて、渦相関法の補正法の誘導について説明する。
【００４８】
以下では、シンチロメータ１によるフラックスおよび消散率にはSASの添え字を、ECに相当する値にはECの添え字を添えて区別することとする。ここで、数式１６〜数式１８において各フラックスの比をとると数式４７〜数式４９となり、シンチロメータ１と渦相関センサによる消散率の比（εｉSAS／εｉEC）と渦相関法によるフラックスより、シンチロメータ１の空間スケールに相当するフラックスが得られる。なお、εｉにおける添字ｉは、何の消散率を表すかによってＴ，ｑ，ｃのいずれかの添字が該当することを意味している。
【数４７】
【数４８】
【数４９】
【００４９】
ここで、εＴSAS／εＴECは、乱流強度（水平風速で標準化した垂直風偏差：σｗ／ｕ）および風向に応答して変化することが分かっている。また、シンチロメータ１のパスに対する風向の変化はシンチロメータ１の空間スケールに影響する。つまり、風向の変化、すなわちパスに対する風の迎え角が変化することで、空間スケール（測定する範囲）が変化する。ただし、このほかに風速や大気の安定度なども影響するが、空間スケールの変化に最も大きく寄与するのが風向と考えられる。そこで、εｉSAS／εｉEC を乱流強度とシンチロメータ１の空間スケールを変数とする関数で表すことで、渦相関センサによる計測値からシンチロメータ１の空間スケールに対応するフラックスを算出することが可能となる。
【００５０】
シンチロメータ１の空間スケールに対応するフラックス算出に関しては以下のとおりである。まず、フラックス測定センサによって検出される乱流信号は風上の領域を起源とし、その広さ（ソースエリア）は風速や大気安定度、そしてセンサの測定パス長に対して変化することになる。シンチロメータ１と渦相関の空間スケールの違いを定量的に解析するために、Kormann and Meixner(2001) のフットプリントモデルを応用してシンチロメータ１のソースエリアを算出すると以下のとおりとなる。すなわち、フットプリントＦ（ｘ）は風上方向にｘ（ｍ）離れた位置のフラックスへの寄与を表した関数で、風速ｕ（ｍｓ−１）、摩擦速度ｕ＊（ｍｓ−１）測定高さｚ（ｍ）、大気安定度パラメータζを変数とする。さらに、De Bruin et al.(2002) の計算例に従い、シンチロメータ１のソースエリア（ＡSAS）を、Ｆ（ｘ）とセンサのパス長で表される横風方向ｙ（ｍ）の感度関数Ｄ（ｘ，ｙ）との積が描く等値線に囲まれた面積から求める（図６参照）。ここで、Ｄ（ｘ，ｙ）はシンチロメータ１が信号を検出する感度がパスの中央部で最も高いことを表す関数で、本実施形態では式形の簡単さから、Wang et al.(1978) が示した感度関数Ｄ（ｘ，ｙ）
［数５０］
Ｄ（ｘ，ｙ）＝（ｙ／Ｌｐ（１−ｙ／Ｌｐ））５／６
を標準化して用いる。ただし、Ｌｐはシンチロメータ１のパス長である。また、風向によって変化するシンチロメータ１のソースエリアは、パスに直交する方角からの風向角（θ）を用いて次のように表すこととする。
【数５１】
シンチロメータ１および渦相関センサのソースエリアの比ＲＡ（＝ＡSAS／ＡEC ）は、風向がシンチロメータ１のパスに平行な場合（θ＝９０°、２７０°）に１となり、直交風向時に増加する（図７参照）。
【００５１】
次に、領域平均フラックスを算出する。ここでは、まず上記で求めたソースエリア比ＲＡおよび乱流強度σｗ／ｕを変数とする関数ｆ（σｗ／ｕ，ＲＡ）の形を定め、一定期間に得られたεＴSAS／εＴEC の実測値に対する回帰によって、関数の係数を決定する。
【数５２】
次に、妥当と考えられるフラックス評価値を導くＲＡを求め、その関係が全ての風向に対して成立すると仮定して、数式４７より顕熱フラックスＨSASを求める。こうして求めたフラックスを領域平均フラックスと定義する。領域平均フラックスは全風向に対してシンチロメータ１の空間スケールを持ち、渦相関センサによる計測値のみでフラックス評価を行うことができるため、シンチロメータ１が測定不能な気象条件のフラックスも評価できる。さらに、水蒸気とＣＯ２が熱と同じ原理で輸送されると考えられることから、
【数５３】
と仮定すると、潜熱とＣＯ２についても同様に領域平均フラックスが計算できる（図８参照）。
【００５２】
なお、上述の実施形態は本発明の好適な実施の一例ではあるがこれに限定されるものではなく、本発明の要旨を逸脱しない範囲において種々変形実施可能である。
【実施例１】
【００５３】
中部地方の浅間東麓の５４年生の落葉広葉樹林においてＣＯ２を含むフラックス測定を、2002年から2004年を通して行った。ここでは、シンチロメータ１の空間スケールから決定した関数を用いて渦相関法の測定値を補正することにより、領域平均フラックスを算出するという「領域平均フラックスの算出」を実施し、算出した領域平均フラックスによって熱収支が正しく評価されていることを確認した。また、光合成量の計測結果とＣＯ２フラックスから求めた群落光合成量とを比較した。優先種であるダケカンバ個葉の光合成特性と葉の分布をもとに算出した森林の光合成ポテンシャルに対し、渦相関法が約ＣＯ２程度過小評価する一方で、領域平均フラックスによる算出値は良好に一致したことから、ＣＯ２交換量においても領域平均フラックスによって評価値の信頼性が向上することが分かった。さらに、森林のＣＯ２吸収量の評価を実施した結果、気象条件や樹種が類似した場所で得られた既往の渦相関法による評価結果とほぼ同等であった。以上につき、以下に実施例１として説明する。
【００５４】
ここでは、高さ２８ｍの鋼製足場材組上げ式のタワーを東西に８６ｍ離して配置し、東側の塔頂で渦相関法を実施するとともに領域平均フラックス測定のために、２基のタワー間をシンチロメータ１の測定パスとした。渦相関法では顕熱、潜熱、ＣＯ２フラックスを測定した。フラックス算出には１０Ｈｚサンプリング時系列の３０分ブロック平均を用いた。ＮＥＥの算出には降雨などの影響を受けにくく、測定値のドリフトが少ない閉光路（closed-path）式ＣＯ２分析計（LI-7000, Licor USA）によって計測したＣＯ２密度変動を用いた。内径４ｍｍのデカボンチューブ４０ｍを経由してタワー基部の観測舎に塔頂の大気サンプルを9.0Ｌmin−１で吸引し、その一部を2.0Ｌmin−１でＣＯ２分析計に導入した。サンプルエアが塔頂から分析器までに到達する時間遅れは、９秒程度であった。サンプルエアーの吸引によって信号の高周波成分が劣化してしまう問題は、乱流のスケールが大きな森林上ではほとんど考慮しなくて良いことが報告されている。当試験林においても、応答の速い開光路（open9-path）式のＣＯ２分析計と比較して、closed-path式のシステムが問題ない精度でＣＯ２フラックスを計測することを確認した。領域平均フラックスを計算するための関数ｆ（σｗ／ｕ，ＲＡ）を定めるために、シンチロメータ１の測定値が得られた2002年６月から１１月までの晴天日の測定値を用いて回帰計算を行った。
【００５５】
次に、欠測データの補間について説明する。フラックスからＮＥＥ（Net Ecosystem Exchange：生態系純交換量）を算出するためには、計測器の故障、メインテナンス、停電などさまざまな理由により発生する欠測を補い、連続した時系列を得る必要がある。そのために、1.5時間以内の短期欠測は線形補完を行い、それ以外はＣＯ２フラックスと相関の高い気象要素を用い、あらかじめ求めておいた気温‐呼吸量および光‐光合成の関係から補間を行った（図９参照）。測定値が得られない場合はAMeDASデータから当試験林の日平均日射量および気温を推定した。
【００５６】
また、フラックス測定値との比較を行うために、葉面積、個葉の光合成量の測定を行った。タワーを中心とした６０ｍ（Ｎ−Ｓ）×１４０ｍ（Ｅ−Ｗ）の範囲を１０ｍ×１０ｍの区画に分け、すべての生態データ取得の基準とした。ここでは、葉面積の季節変動を明らかにするために、開口部１ｍ×１ｍの落葉採取装置（リタートラップ）を林床から1.3ｍの高さに千鳥状に１９個配置し、９月中旬から１１月上旬の落葉期に落葉・落枝を定期的に採取した。また、強風などによる落葉期以外の落葉・落枝を考慮して着葉期の６月から９月上旬まで、平均的な落葉量を示した３個所のリタートラップを用いて落葉量を調査した。あらかじめ求めておいた落葉の乾燥重量と面積の１次回帰式から、採取した落葉サンプルの葉面積を求め、全葉面積から調査区域のＬＡＩ（Leaf Area Index：波面積指数）を求めた。光合成量の測定は、渦相関法を実施したタワーの北側に隣接したダケカンバ（Betula ermanii Cham.）の高度９ｍ付近にある任意の個葉を対象に、携帯型光合成蒸散量測定装置（LI-6400, Licor USA）を用いて行った。ダケカンバは、事前の毎木調査によって、当試験林を構成する樹木の地上部（枝＋幹）炭素量の44.3％を占める優先種であることが分かっている。当装置は小型のリーフチャンバを備え、チャンバ内の気温、湿度およびＣＯ２濃度を外部環境に追従させて制御できる。ダケカンバ葉が十分に光飽和する光強度から７段階（2000, 1500, 1000, 500, 200, 100, 0μmolm-2s-1）に変化させ、最も光合成能力が高いと思われる７月から９月期間で光強度に対する光合成量の反応（光−光合成反応）を調べた。
【００５７】
以上の結果、以下のような結果が得られた。まず、2002年から2004年にかけての年毎のフラックスデータ欠測率はそれぞれ１８％、５％、３％であった。2002年に欠測率が高いのは落雷による停電の影響で、そのほかにはメインテナンスや悪天候、機器の故障などの原因があった。ここで、FLUXNET(Falge et al.,2002)や、Asia Flux(Asia Flux運営委員会、2003) などでは信頼性の高いフラックス測定のために、乱流データのクオリティコントロール（ＱＣ）が提言されている。ＱＣでは乱流の成立性が乏しい状況での渦相関法は信頼性が低いとして、摩擦速度ｕ＊や時系列データの定常性に基準を設け、基準に満たないデータを棄却して補間の対象とする。このため、標準的な森林フラックスサイトの欠測率は３０％程度と言われている。しかし、本実施例では、渦相関法の信頼性が低いとされる乱流の乏しい状態あるいは対流の大きな状態における計測値も含めてシンチロメータ１のデータと比較し、その結果に基づいて領域平均フラックスを算出するため、ＱＣは適用せずすべてのデータを用いた。
【００５８】
続いて、補正関数の決定を行った。ここでは、シンチロメータ１と渦相関センサによる温度変動の消散率比εＴSAS／εＴEC は、機械的な乱流強度（σｗ／ｕ）が大きな中立時では、ＥＣセンサによる消散率がＳＡＳを上回り、（σｗ／ｕ）が小さな安定時および不安定時ではＳＡＳの消散率がＥＣより大きくなった（図１０（ａ）、図１０（ｂ）参照）。消散率とフラックスの絶対値は数式４７のように比例関係にあるため、この結果は、対流が卓越する条件や成層した大気条件では渦相関センサがフラックスを過小評価するというKanda et al.(2004)の結果と一致した。また、シンチロメータ１のソースエリアが渦相関センサに比べて大きくなる（ＲＡが大きくなる）につれてεＴSAS／εＴEC の分散の上端が大きくなる傾向が見られた（図１０（ｃ）、図１０（ｄ）参照）。センサの空間スケールが大きくなることで、より広い領域を起源とする信号の平均値を検出できるため、空間的スケールの大きな大気の乱流現象を捉えることができ、結果的に上記のσｗ／ｕに対するεＴSAS／εＴEC の関係が強調されると考えられた。
【００５９】
そこで、図１０（ａ）および図１０（ｂ）のｆ（σｗ／ｕ，ＲＡ）が大きい領域の漸近値をオフセットと考え、指数減衰関数（数式５４、数式５５）によって上記の関係を表現することとした（図１０の右欄参照）。ここでは、数式中の係数ａ１，ａ２，ａ３を昼間のおよび夜間について最小二乗法によって決定した。夜間と昼間で係数が異なるのは、ソースエリアの算出に用いたフットプリントモデルが夜間と昼間で異なる関数で表されるためである。
［数５４］
ｆ（σｗ／ｕ，ＲＡ）＝ｙ０＋ａ１exp［−（σｗ／ｕ−ｘ０）／ａ２］
［数５５］
ｘ０＝（ＲＡ−１）／ａ３
【００６０】
関数ｆ（σｗ／ｕ，ＲＡ）と消散率比εＴSAS／εＴEC の比較を図１２に示す。昼間（図１２（ａ））に比べて夜間（図１２（ｂ））の分散が大きいのは、渦相関センサによる温度変動の消散率εＴEC の算出精度が乱流状態の乏しい状況で低下したことが原因と考えられた。
【００６１】
次に、領域平均フラックスが得られると考えられる最適なＲＡを計算するために、年間の熱フラックスの収支とＲＡの関係を調べた（図１３）。昼間の不安定時および夜間安定時共に、ＲＡの増加にしたがってインバランス率Ｉ
［数５６］
Ｉ＝（Ｈ＋λＥ−Ｒｎ）／Ｒｎ
が０に近づき、熱収支が釣り合う効果が見られた。特に、夜間安定時はその効果が強かった。潜熱と顕熱の他に林冠内の大気や植物体に蓄えられる熱や地中伝熱（ 林冠内貯熱）が熱収支に関与するが、その値は大きくても１０数Ｗｍ−２であることと、ＲＡが２０を超える測定値がほとんど無く、外挿による誤差の拡大を防ぐ観点から、インバランス率＝−0.3を示すＲＡを求めた（夜間10.8、昼間23.1）。これによってｆ（σｗ／ｕ，ＲＡ）がσｗ／ｕの変数として表され、渦相関センサの測定値に基づいて領域平均フラックスが得られた。
【００６２】
ここで、領域平均フラックスの特徴について記す。領域平均フラックスは、σｗ／ｕが小さな領域では渦相関法によるフラックスより大きな値を示し、σｗ／ｕが大きな領域では渦相関法より小さな値を示す。その結果、夜間の放出および昼間の吸収量ピークが大きく見積もられ、日変動の振幅が拡大する。図１４に熱およびＣＯ２フラックスの測定例を示す。熱フラックスでは領域平均フラックスによる熱収支の改善が見られる一方で、大気が中立な条件では、領域平均フラックスが渦相関法を下回る状況も見られた（図１４のＡ参照）。また、熱収支の改善効果が小さい状況も見られた。林冠内貯熱の寄与が大きい時間帯が存在する可能性もあると考えられた。ＦＣについても同様に放出・吸収ピークが拡大する傾向が見られ（図１４のＢ参照）、ＦＣに基づいて計算される群落光合成量や、群落呼吸量も変化することになる。
【００６３】
さらに、領域平均フラックスと土壌呼吸量との比較を行った。土壌中の動物および根の呼吸と微生物による有機物の分解によってＣＯ２が放出される土壌呼吸と、地上部の植物体による呼吸が合算して森林からのＣＯ２放出となり、その放出量を群落呼吸量と呼ぶ。森林上で観測したＣＯ２からは、群落呼吸量に相当する情報が得られる。一方、土壌呼吸量は、チャンバ法などによって実測が可能である。そこで、ＣＯ２フラックスの妥当性を検証するために、群落呼吸量と土壌呼吸量を比較した。
【００６４】
まず、群落光合成量Ａｇ、群落呼吸量ＲｅとＮＥＥの間には
［数５７］
ＮＥＥ＝−Ａｇ＋Ｒｅ
の関係がある。夜間はＡｇ＝０なので、ＮＥＥが群落呼吸量となる。群落呼吸量は温度の相関が高く、温度に対する群落呼吸量のパラメタリゼーションには数式５８が用いられる（Lloyd and Taylor,1969）。
【数５８】
【００６５】
ここで、Ｔ(Ｋ)は温度、Ｔ０(Ｋ)は任意に決定する基準温度で、Ｒ０はＴ０(Ｋ)におけるＣＯ２フラックスを表す。また、Ｑ１０は温度係数で、基準温度におけるＣＯ２フラックスに対する、基準温度から１０Ｋ上昇したときのＣＯ２フラックスの比を表す。気温と群落呼吸量の関係（2004年）によると（図１５参照）、見易さのため σｗ／ｕ＞0.3以上のデータを示した。そこで、当該森林の群落呼吸量を解析するために、ＳＡＳおよびＥＣから求められる夜間のＣＯ２フラックス（ＦＣ）のデータと温度の関係から、パラメータＱ１０，Ｒ０を決定した（表１参照）。ここでは、冬季間に積雪することを考慮して、Ｔには地温ではなく地上高さ３ｍの位置の気温を用いた。なお、群落呼吸量Ｒｅは、土壌呼吸量（土壌中の生物による呼吸と根の呼吸）および植物地上部の呼吸量を含む。
【表１】
【００６６】
渦相関法にもとづく群落呼吸量のパラメータＱ１０と領域平均フラックスによるＱ１０は年間の差異において異なる傾向を示したが、いずれも平均気温との間には明確な傾向は見られなかった。次に、夜間の温度−呼吸量の関係が昼間でも同様に成り立つと仮定して、昼間の気温と数式５８から群落呼吸量を推定した。こうして求めた温度と群落呼吸量の関係と後述する群落光合成量とを用いて、数式５７から欠測の補間を行った。以前、土壌呼吸チャンバを用いて評価した当試験林の2002年における年間の土壌呼吸量は6.6tCha−１yr−１であったのに対し、渦相関法と領域平均フラックスにもとづく群落呼吸量はそれぞれ、7.4tCha−１yr−１と7.3tCha−１yr−１で（表２参照）、フラックス測定手法間に明確な相違は見られなかった。一方、2003年および2004年では、渦相関法の3.9tCha−１yr−１および7.2tCha−１yr−１に対して、領域平均フラックスは6.7tCha−１yr−１および8.9tCha−１yr−１と大きな値を示した。群落呼吸量は、土壌呼吸量に植物地上部の呼吸量が加わるため、フラックスによる評価値がチャンバー法を上回るのは妥当な結果といえる。ただし、温度に対する群落呼吸量の分散が大きく（図１５参照）、カーブフィッティングの際の誤差が大きいこととなった。領域平均フラックスから算出した群落呼吸量は土壌呼吸チャンバの測定結果と一致し、領域平均フラックスによって正しく群落呼吸量が評価されることが確かめられた。
【表２】
【００６７】
続いて、領域平均フラックスと群落光合成量との比較を行った。本実施例では、ＣＯ２フラックスの妥当性を検証するために、群落光合成量との比較を行った。ＣＯ２フラックスによって評価される群落光合成量は、森林全体の光合成量を表す。一方、葉の光合成量はリーフチャンバによって実測が可能である。しかし、森林は複雑な立体構造を持ち、林床面積あたりにその数倍の面積の葉が存在する。このため、個葉において計測された光合成量を群落光合成量と比較するためには個葉の計測結果を群落に拡張する必要がある。群落光合成量は、群落最上端に位置する葉の最大光合成量と葉面積から評価することができる。
【００６８】
まず第一に、群落光合成量のラメタリゼーションを行った。ＮＥＥと群落呼吸量の差は同化作用による群落光合成量（数式５９中のＡｇ）を示す。群落光合成量（μmolm-2s-1）Ｐｍは最大光合成量（μmolm-2s-1）を表し、αは初期光合成効率を表す。光合成の単位にはＣＯ２の分子量を乗じた密度表現（mgＣＯ２m-2s-1）も用いられる。数式５７に従って入射光量と数式５８に昼間の気温を代入して求めた呼吸量を差し引いたＦＣのプロットに数式５９をフィットしてパラメータを決定した。こうして求めたパラメータと入射光量および気温からＡｇとＲｅを算出し、欠測の補間を行った。
【数５９】
【００６９】
次に、個葉の光合成量を群落光合成量Ｐｃａｎに換算するための葉面積の算出結果を示す。落葉量から積算した葉面積指数（ＬＡＩ（Leaf Area Index：波面積指数））は、2002年から2004年までそれぞれ5.6，6.3，6.8で、増加傾向を示した。
【００７０】
また、図１６に示した葉面積指数の変化を見ると、展葉が完了した６月から、落葉が開始する９月までの期間も落葉が継続している。当試験林の優先種であるダケカンバは、春の展葉の後に２度目の伸張をすることが認められた。６月から９月の期間中も２次伸張による展葉が継続している場合、その間の落葉量をＬＡＩに算入するとＬＡＩを過大評価することになる。林床面積あたりに存在する実際の葉面積に近い値として、落葉開始直前の９月上旬の値を用いると、それぞれ、5.4，5.6，5.8となり、年間の格差は小さくなったが大きさの順位は変わらなかった。全落葉重量に占めるダケカンバの落葉重量は38.2％，42.1％，42.7％とほぼ41％程度であり、地上部炭素量の割合44.3％と近い値を示した。葉面積換算では29.1％，31.9％，32.5％とほぼ31％程度であった。
【００７１】
続いて、群落光合成量の算出を行った。領域平均フラックスによるＰｍは、７月および８月に最大値を示し、その値は1.7，1.9，1.8（mgＣＯ２m-2s-1）であった。一方、リーフチャンバによって計測したダケカンバ個葉の最大光合成量の平均値はそれぞれ0.8±0.2，0.9±0.1，0.9±0.1（mgＣＯ２m-2s-1）であった。リーフチャンバで得られたダケカンバ個葉の最大光合成量を群落最上端の最大光合成量とみなすと、Ｂｅｅｒの法則に従い、積算葉面積Ｆの位置で葉群に吸収された光量Ｉｉｎｋは以下で表される
【数６０】
ここで、Ｉ０は群落上端の入射光、Ｋは葉の角度などで決まる群落内の吸光定数で、0.5とする。ａは葉の透過率’（＝0.14）である。次に、入射光Ｉｉｎｋに対する光合成量Ｐは数式５９のPARをＩｉｎｋとすることで算出できる。図１７に、積算葉面積に対する葉群の吸光量Ｉｉｎｋおよび数式５９にＩｉｎｋを代入して求めた光合成量Ｐを示す。そして、ＰをＦで積分することにより、群落全体の光合成量Ｐｃａｎを計算することができる。
【数６１】
群落下層の植物は群落上端に比べて光合成能力が低い（最大光合成量が小さい）のが一般的で、実際にはさらに群落光合成量が小さくなることが考えられた。また、このモデルでは水分不足による気孔閉塞がもたらす昼寝現象などは考慮していない。したがって、本計算結果は群落光合成のポテンシャルと考えるのが適当である。植物体内の水分状態が気孔の開閉に影響したり、乱流の状態によってＣＯ２の交換量が変化するなど、フラックスの瞬時値はほぼ群落光合成のポテンシャル以下に分布すると考えられる。図１８に示すように、領域平均フラックスによるＡｇの分布の上端とＰｃａｎがほぼ重なる。このことは、領域平均フラックスが当該森林の光合成特性をよく表しているといえる。
【００７２】
また、図１９に当試験林で観測される最大光合成有効放射量2000μmolm-2s-1における群落光合成のポテンシャルとフラックスから導いた群落光合成量を示した。渦相関法の過小評価はＣＯ２フラックスにおいても明らかで、領域平均フラックスによる最大光合成量は、ＥＣによる最大光合成量の約1.5倍と、理論値である群落光合成のポテンシャルに近づいた。この結果、領域平均フラックスによって森林の光合成が正しく評価されることが示された。
【００７３】
続いて、領域平均フラックスの妥当性評価を行った。本発明者は、以前、シンチロメータ１の適用による熱収支式の整合をフラックス測定の妥当性の根拠とし、同じ原理で測定されるＣＯ２についても、フラックスが妥当に測定されると判断したが、本実施例においては更に領域平均フラックスから導いた群落呼吸量および群落光合成量と、群落呼吸量および個葉の光合成特性に基づく群落光合成量を比較した結果、領域平均フラックスによる算出値は、土壌呼吸量および光合成量の実測値とよく一致した。これにより、領域平均フラックスによるＣＯ２評価の妥当性が示された。また、熱フラックスの収支を判断材料として、フラックス測定の妥当性を評価する方法論が支持される結果となった。以上のような本実施例における結果は、当試験地特有の現象を捉えたものであるが、既存のフラックス測定サイトおよび新規測定サイトにおいても、渦相関法とシンチロメータ１と組み合わせた観測によってサイト固有の補正関数を求め、領域平均フラックスを求めることが可能と考えられた。
【００７４】
さらに、領域平均フラックスに基づく森林ＣＯ２吸収量の評価を行った。表２に示したように、ＮＥＥは大きな年変動を示した。これには、降水量、日射量、気温および水蒸気の飽差などの気象条件に加えて、窒素負荷や病害虫なども複合で影響すると考えられた。例えば、2004年は他の年より落葉時期が早く、９月中に葉量を減らしたにもかかわらず生産量は高かった。領域平均フラックスによるＮＥＥは渦相関法を上回り（図２０のＡ参照）、その結果、群落呼吸量および群落最大光合成量も大きく見積もられた（図２０のＢ参照）。３年間のＮＥＥの平均値は渦相関法では180gCm-2yr-1、領域平均フラックスでは266gCm-2yr-1と約1.5倍になった。以前、当試験林と類似した立地条件の４０年生のダケカンバ林で1998年から1999年にかけて実施された渦相関法によるフラックス測定では163gCm-2yr-1で、本実施例の結果に近い値であった。また、熱収支のインバランスを解消させるようにフラックスを補正したところ、214gCm-2yr-1となり約1.3倍となった。生態系呼吸量（Ｒｅ）と生態系純生産量（ＮＥＥ）の和で表されるＧＰＰ(Gross primary production)について比較すると、従来法で875gCm-2yr-1、補正後では1146gCm-2yr-1で、当試験林の３年間平均値811gCm-2yr-1および領域平均フラックスによる1156gCm-2yr-1と同等の結果であった。これは本実施例とは異なった補正法を用いているが、本実施例の結果と同様に渦相関法による過小評価が起こっていると考えられた。
【００７５】
以上のとおり、本実施例においては、森林におけるＣＯ２吸収量を高い精度で評価することを目的として、シンチロメータ１の空間スケールから決定した関数を用いて渦相関法の測定値を補正することにより、領域平均フラックスを算出した。ここでは、熱収支が整合するように調整した領域平均フラックスと、森林における群落呼吸量や群落光合成量との比較を通して領域平均フラックスの確かさを検証した。その結果、領域平均フラックスを基に算出した群落呼吸量は、以前に開発されたチャンバ法による土壌呼吸量の評価結果との一致を見た。また、従来法である渦相関法から算出した群落光合成量がダケカンバ個葉の計測結果から導いた群落光合成量を大きく下回るのに対し、領域平均フラックスによる群落光合成量は個葉の計測結果とよく対応することが明らかになった。以上より、森林における熱収支および群落のＣＯ２収支において、本発明者が領域平均フラックス評価法として開発した手法が従来の渦相関法の信頼性を高めることが分かった。例えば、従来法によって評価したＣＯ２吸収量の３年間の平均値は180gCm-2yr-1で、類似の立地条件にある落葉広葉樹林の値と同等であった。それに対し、本発明にかかる測定手法を用いて推定した森林によるＣＯ２吸収量の３年間の平均値は266gCm-2yr-1で、渦相関法による評価値の1.5倍であった。
【実施例２】
【００７６】
風速、乱流強度、ＴＫＥ（乱流運動エネルギー）、温度、比湿、ＣＯ２の各要素に関し、シンチロメータ１（ＳＣとも表記している）および渦相関センサ（ＳＣ）間の消散率比と風向の関係について実験を行った（図２１参照）。図中の角度は方位を表し、北を０°としている。グラフ中の値が大きければ大きいほどシンチロメータ１による消散率の比が大きく、例えば風速でいえば、北風ほどシンチロメータ１による消散率の比が大きくなる傾向が強い、という結果が得られた。
【実施例３】
【００７７】
摩擦速度ｕ＊、顕熱フラックスＨ、潜熱のフラックスλＥ、そしてＣＯ２フラックスｆｃのそれぞれに関し、シンチロメータ１および渦相関センサによるフラックス算出値の比較を行った（図２２参照）。図中、四角い記号で示すシンチロメータ１のパスに平行な西風、丸い記号でシンチロメータ１のパスに直交する北風のデータを示した。各期につき、シンチロメータ１による測定結果と渦相関センサによる測定結果の比は図に示すように、シンチロメータ１による測定結果が大きくなる傾向にあった。
【実施例４】
【００７８】
シンチロメータ１と渦相関センサのそれぞれの場合の非着葉期におけるインバランス率をグラフに表示することにより、熱収支インバランスが改善されたことの確認を行った（図２３参照）。インバランス率ＲＩ、熱収支式はそれぞれ図中に示す式で求めることができる。この結果、本発明にかかるＣＯ２フラックス測定法によれば従来法よりもインバランス率が改善されることがわかった（図２３参照）。
【図面の簡単な説明】
【００７９】
【図１】本発明のフラックス測定法におけるフラックスの定義と測定法を説明するための概略図である。
【図２】シンチロメータの測定原理を示す図である。
【図３】本発明におけるフラックス算出のデータ処理の流れを示す図である。
【図４】シンチロメータのパス長による測定領域の違いを渦相関センサの場合も併せて説明するための図である。
【図５】森林生態系のＣＯ２交換の概念を示す図である。
【図６】風向がパスに直交するときのフットプリントおよびソースエリアを示す図である。
【図７】ソースエリア比ＲＡと風向θの関係を示す図である。
【図８】領域平均フラックス算出の手順を示す図である。
【図９】欠測データの補間のフローを示す図である。
【図１０】εＴSAS／εＴEC に対する乱流強度σｗ／ｕおよびソースエリア比ＲＡの関係、ならびに補正関係ｆ（σｗ／ｕ，ＲＡ）と乱流強度σｗ／ｕおよびソースエリア比ＲＡの関係を示す図である。
【図１１】構造関数の例と曲線フィットによる消散率決定の内容（2002年６月10日）を示す図である。
【図１２】εＴSAS／εＴEC とｆ（σｗ／ｕ，ＲＡ）の比較を示す図で、ａ）は昼間、ｂ）は夜間を表している。
【図１３】ソースエリア比ＲＡとインバランス率Ｉ（ただし、Ｉは３年間の平均値）の関係を示す図である。
【図１４】熱フラックスおよびＣＯ２フラックスの日変化における渦相関法ＥＣと領域平均フラックスＳＡＳの関係を示す図である。
【図１５】気温と群落呼吸量の関係（2004年）を示す図であって、見やすさのためσｗ／ｕ＞0.3 の条件部分についてのみ表しているものである。
【図１６】リター積算に基づいた落葉期のＬＡＩの変化を示す図である。
【図１７】葉群の吸光量Ｉ０および光合成量（Ｐ）の積算葉面積（Ｆ）に対するプロファイルを示す図である。
【図１８】群落光合成のポテンシャルおよび群落光合成量の入射光量（PAR）との関係（2004年７月）を示す図である。
【図１９】群落光合成ポテンシャルと群落光合成量の比較（PAR＝2000μmolm-2s-1）を示す図である。
【図２０】領域平均フラックスと渦相関フラックスによる生態系純交換量（ＮＥＥ）および群落最大光合成量（Ａｇ）の比較を示す図である。
【図２１】シンチロメータ（ＳＣ）および渦相関センサ（ＥＣ）間の消散率比と風向の関係を各要素について示す図である。
【図２２】摩擦速度ｕ＊、顕熱フラックスＨ、潜熱のフラックスλＥ、ＣＯ２フラックスｆｃのそれぞれに関し、シンチロメータおよび渦相関センサによるフラックス算出値の比較結果を示す図である。
【図２３】シンチロメータと渦相関センサのそれぞれの場合の非着葉期におけるインバランス率が改善されたことを示す図である。
【符号の説明】
【００８０】
１ シンチロメータ
２ トランスミッタ
３ レシーバ
【技術分野】
【０００１】
本発明は、森林のＣＯ２フラックス測定法に関する。さらに詳述すると、本発明は、渦相関法（ＥＣ）により評価されていた森林全体のＣＯ２吸収量の評価技術の改良に関する。
【背景技術】
【０００２】
温暖化ガス削減策として京都メカニズムにおける森林吸収源を利用するには、森林吸収源におけるＣＯ２の吸収量の評価方法を示し、実行する必要がある。温暖化抑制や気候変動における森林の影響を定量化する上では、森林全体のＣＯ２収支の評価（フルカーボンアカウンティング）が重要である。このような森林全体のＣＯ２収支は、ＣＯ２の交換量（フラックス）を積算することで評価できるものであり、現在、地球規模のＣＯ２収支の評価・予測が盛んに行われている。このように評価や予測を行うにあたっては、地表面で実測されたエネルギー・物質交換量（フラックス）の検証データとしての利用価値が高いことから、フラックスサイトをネットワーク化し、広域的かつ長期間にわたってフラックスデータをグラウンドトゥルースとして利用できる環境の整備が世界的に進められている。また、現在、ＩＰＣＣ（Intergovernmental Panel on Climatic Change、1988年に設立された国連の組織で、地球温暖化問題についての科学的な議論を行う場）が推奨するＣＯ２吸収量の評価手法があるが、この評価手法は森林のＣＯ２収支のうち土壌からのＣＯ２放出を正確に把握できないことから、森林のＣＯ２収支を正しく評価できる新しい手法が今後は求められるであろうという面もある。
【０００３】
例えば、小規模な公園程度の大きさであればＣＯ２の出入りの量に差が生じるが、数十ｈａあるいは数百ｈａといった大規模な森林であれば出入りの量はほぼ等しいと考えることができる。したがって、森林におけるＣＯ２吸収量を評価するにあたり、森林が横方向に均一だと考えれば、森林に対して横方向から流入する風と横方向に流出していく風のＣＯ２の流量（交換量）はほぼ等しいとして扱うことが可能となる（図１参照）。したがって、この考えに基づけば、横方向の流量は考慮せず、森林の上部空間において上下方向に出入りする空気量などを例えば１年間測定することにより、１年の間に森林がどの程度ＣＯ２を吸収したか、あるいは放出したかがわかる（図１参照）。より具体的には、森林の上下方向におけるＣＯ２の交換量速度を把握し続け、最後に積算すればある期間（例えば１年）の吸収量（放出量）がわかる。この吸収ないしは放出の量、つまりはＣＯ２の交換量を単位時間あたりの量として考えたものが「フラックス」、つまりいわば交換速度ということになる（図１参照）。
【０００４】
ここで、森林上におけるＣＯ２の収支を表すＮＥＥ（Net Ecosystem Exchange：生態系純交換量）は、ＣＯ２フラックスＦＣを用いて次式で表される。
［数７］
ＮＥＥ＝ＦＣ＋ＦＳ
ＦＳは森林内の空間に一時的に存在するＣＯ２量（林冠内貯留量）で、１日の中で貯留と放出を繰り返し、１日を通した収支は非常に小さい。したがって、１年間のＮＥＥは、ＣＯ２フラックスの時系列１年分を積算することで求めることができる。１年間のＮＥＥ（mgＣＯ２m-2yr-1）は光合成の同化作用によって森林生態系が１年間に吸収したＣＯ２量を示す（図５参照）。
【０００５】
従来、このＣＯ２フラックスを測定するための手法の一つとして、渦相関法（ＥＣ：Eddy-covariance method）が広く用いられている。この渦相関法は、垂直風速の変動（ｗ’）と気温や湿度、ＣＯ２密度変動の測定値（ｃ’）の共分散の平均値（ｗ’ｃ’（文中ではアッパーラインを付していないが、算術平均を表しているものとする））から乱流輸送量であるフラックスを求める、つまり、図１に示してあるように、
ＣＯ２フラックス＝風の乱れ×ＣＯ２密度変動＝ｗ’ｃ’
という式に従い、ＣＯ２フラックスを求めるというものである。この渦相関法による測定手法は、各物理量を直接測定するものであることから最も仮定が少ない測定法とされている。また、この渦相関法を利用する場合の施設としては森林上空の風の流れとＣＯ２の密度の変動を直接的に測定するための施設、例えば森林内に建てられたタワーや当該タワーに設けられたセンサなどが利用されている（例えば、非特許文献１参照）。センサとしては、風の乱れを測定するための超音波風速計や、ＣＯ２密度変動を測定するためのＣＯ２／Ｈ２Ｏ分析計などが用いられる。
【０００６】
【非特許文献１】Thiermann, V. "A displaced-beam scintillometer for line-averaged measurements of surface layer turbulence". Tenth symposium of turbulence and diffusion, 29 Sept -2 Oct, 1992, Portland, OR., published by the American Meteorological Society, Boston, MA.: ｐ244-ｐ247 (1992)
【発明の開示】
【発明が解決しようとする課題】
【０００７】
しかしながら、従来の測定法の場合、測定精度が十分ではないとの見方がある。すなわち、フラックスは大気下層の地表面付近では大気乱流によって駆動されることから、晴天日の夜間など大気が安定した状態や、大気の対流が大きい不安定な条件では、渦相関法による熱フラックスと放射計による放射フラックスが釣り合わない現象（インバランス）が多発している。そして、このような現象（インバランス）は、主に渦相関法に起因する問題と認識されている。
【０００８】
そこで、本発明は、森林全体におけるＣＯ２の吸収量を高精度に評価するための森林のＣＯ２フラックス測定法を提供することを目的とする。
【課題を解決するための手段】
【０００９】
かかる目的を達成するため、本発明者は種々の検討を行った。森林全体におけるＣＯ２の吸収量を高精度に評価するため、森林全体における時間平均ないしは空間平均のＣＯ２フラックスを測定することを考えた。ただし、上述したように、渦相関法だと測定精度が十分ではないとの見方があり、本発明者はこの問題について検討した結果、森林上空の１点にて測定を行ういわば「点測定」型であることに理由があるのではないかと考え、さらに検討を進めた結果、パス（トランスミッタとレシーバとの設置間隔）を長くとれるシンチロメータを使ってＣＯ２フラックスをいわば面的に広い領域にて測定するという着想を得るに至った。シンチロメータならば、トランスミッタとレシーバとを数十ｍあるいは100ｍ以上離して設置し、それらの間におけるレーザ光の屈折率変化を測定することにより、少なくとも渦相関法の場合よりも広いエリアについてＣＯ２フラックスを測定することが可能となる。
【００１０】
ただし、シンチロメータの場合だと、温度の乱流強度しか測定することができないために、上述したＣＯ２フラックスを求めるための要素である「ＣＯ２密度変動」に関するデータを直接求めることができないという欠点もある。つまり、より広い領域にてフラックス測定することにより測定精度を高めることが期待できても、ＣＯ２フラックスを求めるための要素を直接求めることができなければそもそもＣＯ２フラックス測定が不可能となってしまう。
【００１１】
ここで、本発明者は、パスを長くとれるシンチロメータを測定手段として利用しつつＣＯ２フラックスを求めるためさらなる検討を加えた。そこで、測定手法の差異に対するソースエリアの影響を調べるために、ソースエリアおよびシンチロメータと渦相関センサによる消散率を比較することに想到し、さらなる検討の結果、ある知見を得るに至った。
【００１２】
本発明はかかる知見に基づくものであり、まず、森林のＣＯ２フラックス測定法として、森林上空のシンチロメータのトランスミッタから照射したレーザ光をシンチロメータのレシーバで受光する際、森林上における乱流あるいは温度変動に起因するシンチレーションに基づく光強度の変動を計測して温度変動の消散率εＴを求め、さらに、森林上空の渦相関センサによって、温度に関する変動の分散の平方根σＴ、水蒸気に関する変動の分散の平方根σｑ、ＣＯ２に関する変動の分散の平方根σｃを測定し、その後、以下の関係式
【数８】
【数９】
に基づき、水蒸気の消散率εｑ、ＣＯ２の消散率εｃをそれぞれ算出し、スカラー量である各消散率εＴ、εｑ、εｃからフラックスを求めることを特徴とするものを発明した。
【００１３】
例えば森林の上空で空気乱流や熱変動が生じると、光の散乱が生じ、光強度が変動する。そこで、このように相関的に変動した光強度に基づき、乱流運動エネルギーの消散率、あるいは温度変動の消散率を求める（図２参照）。また、空気の乱れとＣＯ２の密度とを掛け合わせる、あるいは空気の乱れと温度変動と掛け合わせることによって熱のフラックスを求めることができる。ここで、フラックスの生成と消散には一対の関係がある。そこで、消散率を求めることによって間接的にフラックスを算出することが可能である。
【００１４】
ただし、この手法においてはＣＯ２のフラックス自体を直接的に求めることはできず、顕熱フラックス（伝熱）と風速シアー（乱流の強さ）という２つが直接的に求められるに過ぎない。そこで、この方法を応用してＣＯ２のフラックスあるいは水蒸気のフラックスを求める手法を新たに導入する。すなわち、上述したようにフラックスの生成と消散とには一対の関係があることから、本発明においては、消散率を求めることによってフラックスを間接的に算出することとしている。この場合、渦相関センサによる計測値からシンチロメータの空間スケールに対応するフラックスを算出する点が特徴的である。
【００１５】
より具体的に説明すると、シンチロメータを用いればフラックスを間接的に測定することは可能だが、このシンチロメータの場合、温度と乱流の強さしかわからない。そこで、温度と水蒸気の消散率の比が、温度と水蒸気の分散の比の二乗と等しいという関係式を導入し、渦相関センサによって計測したデータから、温度に関する変動の分散の平方根σＴと水蒸気に関する変動の分散の平方根σｑの比を求め、さらには温度に関する変動の分散の平方根σＴとＣＯ２に関する変動の分散の平方根σｃの比を求め、ここから、シンチロメータによっては直接的には測りえない水蒸気の消散率εｑ、ＣＯ２の消散率εｃをそれぞれ算出することとしている（図３参照）。こうしてスカラー量である各消散率εＴ、εｑ、εｃを求めることができれば、これらの結果からフラックス（例えばＣＯ２フラックス）を求めることが可能となる。
【００１６】
なお、渦相関法において熱フラックスの過小評価が見られたことがあった。このような測定手法に対し、本発明の場合、シンチロメータによる熱フラックスは熱収支のインバランスを緩和させ、より妥当な評価結果が得られる。また、ＣＯ２フラックスは顕熱および潜熱フラックスと同じ測定原理に基づくことから、シンチロメータの適用によって、より信頼できる値が得られることになる。
【００１７】
ここで、本発明者は検討を重ねた結果、森林上に適用したシンチロメータから得られた乱流運動エネルギーおよびスカラー分散の消散率（フラックス算出のための支配的パラメータ）が、測定パスに直交する風向時に渦相関センサによる消散率を特徴的に上回る現象を知見するとともに、この現象の原因のひとつに、シンチロメータと渦相関センサの間でソースエリアが大きく異なることがあると考えた。かかる知見に基づいて、上記のようなＣＯ２フラックス測定を実施する場合には、シンチロメータにおけるレーザ光の照射方向たるパスの方向が森林における風向と直交するようにトランスミッタおよびレシーバを配置しまたは配置を変換することが好ましい。シンチロメータによる測定領域を広げるという観点からは、上記のようにパス方向を風向に直交させれば、シンチロメータを利用してＣＯ２フラックスを測定する場合により大きな効果を得ることが可能となる（図４参照）。例えば、風向が一定ないしは安定となっている場合にパスを直交させることとすれば測定精度が高くなる。また、特に測定パスをより長距離化した場合に、シンチロメータは広い領域においてフラックスが導出できるため、より広い面的な測定が可能となるという利点がある。本発明の場合、トランスミッタとレシーバの少なくとも一方を移動可能としておけば風向に応じてパスの方向を適宜変化させることが可能となる。
【００１８】
以上がシンチロメータを利用した森林のＣＯ２フラックス測定法であるが、本発明者はさらに検討を重ねた。シンチロメータを利用すれば上記のような効果が得られる反面、レーザ光の光路が遮られると測定不能となるため、雨天や霧など悪天候時のデータが得られないという面もある。また、シンチロメータの測定パスに風向が直交する際には領域を代表する値が得られる一方で、パスに平行な風向時には渦相関法と同等の評価値となる。フラックスの積算によってＣＯ２吸収量を得るためには、シンチロメータの空間スケールを持つフラックスを長期に安定して算出する必要もあると考えた。そこでさらに検討を重ねた結果、シンチロメータと渦相関法による温度変動の消散率の比εＴSAS／εＴECが、(i)センサの測定範囲の比（換言すればソースエリアの比、ＲＡ＝ＡSAS／ＡEC ）および (ii)大気の乱流強度（σｗ／ｕ） と相関が高いことから、これら２つの要素をパラメータとして、渦相関法の測定値をシンチロメータ相当のフラックスに変換する手法に想到するに至った。
【００１９】
請求項１に記載の発明はこのような考えに基づき想到したものであり、森林全体のＣＯ２吸収量を測定するためのＣＯ２フラックス測定法において、シンチロメータを用いた場合と渦相関法による場合の温度変動の消散率の比εＴSAS／εＴECが、シンチロメータを用いた場合の測定範囲ＡSASと渦相関センサによる場合の測定範囲ＡEC の比ＲＡ（＝ＡSAS／ＡEC ）との相関、および大気の乱流強度σｗ／ｕ との相関が高いことに基づき、これら２つの要素ＲＡ（＝ＡSAS／ＡEC ）およびσｗ／ｕ をパラメータとする関係式
【数１０】
を、シンチロメータおよび渦相関センサを使って実際に得られた一定期間における測定結果に基づきあらかじめ用意しておくとともに、数式１０中における測定範囲の比ＲＡには一定値を与えることによりすべての風向に対するシンチロメータの広い測定範囲を適用することとし、その後、森林の上空に実際に設置した渦相関センサによって大気の乱流強度σｗ／ｕを測定し、以下の関係式
【数１１】
に基づき、渦相関法によって求めた顕熱フラックスＨECをシンチロメータによって求めた場合に相当する顕熱フラックスＨSASに補正し、尚かつ、渦相関法によって求めた潜熱フラックスλＥECおよびＣＯ２フラックスＦｃECに関しても、以下の関係式
【数１２】
【数１３】
に基づき、シンチロメータによって求めた場合に相当する潜熱フラックスλＥSASおよびＣＯ２フラックスＦｃSASに補正し、シンチロメータにて測定したのと同様の測定結果を得ることを特徴としている。この場合のＲＡは風向によって変化するシンチロメータの空間スケール、σｗ／ｕは大気の乱流状態を表している。
【００２０】
本発明の場合、上述の(i)と(ii)を変数とする関数（数式１０）を、実際にシンチロメータおよび渦相関センサを使った過去の測定結果に基づきいわば経験値としてあらかじめ用意しておく。その後の計測は渦相関法のみを実施し、シンチロメータを使った実測はもはや行わないが、あらかじめ用意された関係式を利用し、フラックスを補正することによってシンチロメータを使った場合と同様の測定結果を得ることができる。つまりは、当該測定対象森林にシンチロメータを実際に設置する必要はなく、あくまでも従来と同様の手順により渦相関センサのみにより実際の測定を実施して測定データを得ることができる。例えば渦相関センサのみによる測定結果は実際の値よりも過小評価される傾向があるが、本発明においては特定の関係式に基づきこの結果を補正ないしは補間し、あたかもシンチロメータを使って得たのと同様、真の値により近い測定結果を算出することができる。
【発明の効果】
【００２１】
前述のＣＯ２フラックス測定法によれば、測定パスを長距離化できるシンチロメータを利用することにより、広い領域においてフラックスを計測することが可能となる。例えば従来の渦相関法の場合、いわば点で測定するため測定領域が極めて狭いのに対し、本発明にかかるフラックス測定法の場合だと大きい領域を測定できるため、狭小領域における変動には依存しない平均的かつ精度の高い測定データを得ることができる。これにより、森林全体のＣＯ２吸収量を高精度に評価することが可能となる。
【００２２】
また、前述のＣＯ２フラックス測定法においては、風向にパスを直交させるか、あるいはこれに近似した状態とすることによってより広い測定領域を確保することができる（図４参照）。これによれば、森林全体のＣＯ２吸収量をより高精度に評価することが可能となる。
【００２３】
請求項１の発明によると、シンチロメータは使わずに、シンチロメータを使ったのと同様の高精度の測定結果を得ることができる。こうした場合、トランスミッタとレシーバのそれぞれを設置するためのタワーを設けなくても済むため、コストを抑えつつ、ＣＯ２吸収量を高精度に評価することが可能となる。要は、測定領域の異なる計測値から導いた補正を施すことによって、あたかもシンチロメータにて計測を実施したのと同様の結果を得ることができ、渦相関法による熱フラックスのインバランスを解消に向かわせることが可能となる。
【００２４】
しかも、シンチロメータは光学的な計測機器であるために雨や霧の場合に稼働率が下がるといったように天候に左右されやすいが、実際にシンチロメータを使わずに高精度な結果が得られるようにした本発明にかかる測定法によれば天候の影響を受けずに随時測定することが可能になるという利点もある。
【発明を実施するための最良の形態】
【００２５】
以下、本発明の構成を図面に示す実施の形態に基づいて詳細に説明する。
【００２６】
図１〜図８に本発明の実施の一形態を示す。本発明にかかるＣＯ２測定法は、森林全体におけるＣＯ２フラックスを計測し評価するというもので、森林上空のシンチロメータ１のトランスミッタ２から照射したレーザ光をシンチロメータ１のレシーバ３で受光する際、森林上における乱流あるいは温度変動に起因するシンチレーションに基づく光強度の変動を計測して温度変動の消散率εＴを求めるか、あるいは経験値としてあらかじめ用意しておき、さらに、森林上空の渦相関センサによって、温度に関する変動の分散（またはその平方根σＴ）、水蒸気に関する変動の分散（またはその平方根σｑ）、ＣＯ２に関する変動の分散（またはその平方根σｃ）を測定し、その後、以下の関係式
【数１４】
【数１５】
に基づき、水蒸気の消散率εｑ、ＣＯ２の消散率εｃをそれぞれ算出し、スカラー量である各消散率εＴ、εｑ、εｃからフラックスを求めるというものである。
【００２７】
まず、ここで本明細書中における各種用語について説明しておく。まず、一般に「フラックス」といった場合、気流の乱れによって鉛直方向に運ばれる乱れのエネルギー（運動量フラックス）や、熱・物質の単位面積あたりの移動速度（顕熱フラックス、潜熱フラックス、ＣＯ２フラックス）のことを意味する。本明細書では、単位時間当たりのＣＯ２などの交換量のことをフラックスと呼ぶ。「渦相関法」は、超音波風速計とＣＯ２／Ｈ２Ｏ変動計を用いて計測した水平および鉛直方向の風速と物質量の変動成分からフラックスを算定する手法である。「消散率」は、大気の乱れのエネルギーや、熱・物質の密度変動の消散を表す値である。接地境界層の乱流エネルギーについての生成と消散に関する収支式から乱流フラックスが算定される。
【００２８】
本実施形態では、まずシンチロメータ１によるフラックスの算出方法の概要について説明する。次に、その手法を拡張し、シンチロメータ１の空間スケールから決定した関数を用いて渦相関法の測定値を補正するというフラックス算出方法について説明する。
【００２９】
シンチロメータ１（本実施形態においてはSAS：Small Aperture Scintillometerを用いている）は、可視レーザーを射出するトランスミッタ２とレーザー強度を検出するレシーバ３の一対で構成されており、乱流によって引き起こされる大気の屈折率の変動を光強度変化として捉え、顕熱フラックス算出のパラメータである温度変動の消散率（εＴ）を算出する（図２参照）。温度変動の消散率は、渦相関センサによって計測した風速と温度の時系列からも計算することができる。フラックスの収支式にモニンオブコフの相似則を適用することで、消散率は次式によってスカラーフラックスである顕熱Ｈ（Ｗｍ−２）、潜熱λＥ（Ｗｍ−２）およびＣＯ２フラックスＦＣ（mgＣＯ２m−２s−１）と関連付けられる。
【数１６】
【数１７】
【数１８】
ここで、ｕ＊は摩擦速度（ｍｓ−１）、ρは空気密度（ｇｍ−３）、ｋｋａｒはカルマン定数（＝0.4）、ｚは測定高さ（ｍ）、ｃｐは空気の定圧比熱（Ｊｇ−１Ｋ−１）およびλは水の蒸発潜熱（Ｊｇ−１）を表す。εｑ、εｃはそれぞれ、比湿（ｑ ｋｇｋｇ−１）変動およびＣＯ２密度（ｍｇＣＯ２ｍ−３）変動の消散率である。φε(ζ)は無次元化された消散率で大気安定度を表すパラメータ
［数１９］
ζ＝ｚ／Ｌｍｏ（Ｌｍｏ：モニンオブコフ長）
の関数として与えられる。
【００３０】
続いて、シンチロメータ１によって決定したｕ＊およびεＴと、同時に測定したｑおよびＣＯ２密度変動の分散から、εｑ，εｃを求めて潜熱フラックスおよびＣＯ２フラックスを算出する手法の内容について説明する。
【００３１】
シンチロメータ１によるフラックス算定は消散法に基づいて行われる。乱れは生成消散によって成り立っており、渦相関法がスカラー変動や乱流運動エネルギー（気流の乱れのエネルギー：Turbulent Kinetic Energy＝ＴＫＥ）生成の直接測定に相当するのに対し、消散法はスカラー変動やＴＫＥの消散率（エネルギーの消散を表す正の値）からそれらの生成を計算する手法である。消散法は船の動揺など低周波ノイズの影響を受けにくいことから海上でのフラックス測定などに応用されてきた。消散法によるフラックス算出では、まずＴＫＥとスカラーの変動の消散率を求め、次いでフラックスと消散率との収支式からフラックスを算定する。シンチロメータ１はこのうちＴＫＥとスカラー変動の消散率を測定する装置である。
【００３２】
消散法の原理について以下に説明する。水平一様で定常な接地境界層における乱流運動エネルギーＴＫＥ(ｅｔ)およびスカラー変動の収支方程式は以下のように書かれる。
【数２０】
【数２１】
【００３３】
ここでｓは気温Ｔ、比湿ｑ、やＣＯ２などのスカラーを表し、ｐは気圧で
［数２２］
θ’ｖ＝θ（1＋0.61ｑ）
は仮温位（ただしθは温位）、ｅｔはＴＫＥ（ｅｔ＝（ｕ２＋ｖ２＋ｗ２）／２）を表す。ここで、ｕ，ｖ，ｗは平均風向、横風方向、鉛直方向の風速成分である。また、この文中ではアッパーラインを付していないが、数式中のバーは平均値、プライムは平均値からの変動を表す。数式２０の各項はそれぞれ生成項、浮力項、輸送項、圧力項、そして消散率である。また、数式２１の各項は、生成項、輸送項、そして消散率である。両式第１項のｕ’ｗ’およびｗ’ｓ’（ただし本文中ではアッパーラインを省略）が求めるべき運動量およびスカラーのフラックスに相当する。接地境界層においては両式とも輸送項および圧力項を無視できると仮定して、数式２３で定義される摩擦速度ｕ＊、
【数２３】
Monin-Obukov長（Lmo）およびKarman定数（ｋkar 0.4）を用いて、それぞれを
ｕ＊3／ｋkar(ｚ−ｄ0) と ｋkar(ｚ−ｄ0)ｕ＊／（ｗ’ｓ’）２ によって無次元化すると、数式２０、数式２１の各項は測定高さｚ(m)において、
【数２４】
【数２５】
と表される。φεとφεｓはＴＫＥおよびスカラーの無次元消散率を、φｍとφｓはＴＫＥおよびスカラーの無次元した鉛直勾配を表す。ｄ0 はゼロ面変位で、植生地上での風速プロファイルの対数分布を仮定した際の、高さを測る基準面を上方に修正するためのパラメータである。群落高×0.7という値がよく使われる。さらに、ここでＴＫＥおよびスカラー変動の消散率は次のように定義される（Brutsaert，1982）。
【００３４】
【数２６】
【数２７】
水平一様な接地境界層ではこれらの無次元消散率は、大気安定度を示すパラメータζ＝（ｚ−ｄ０）／Ｌｍｏ（安定時ζ＞０−、不安定時ζ＜０）の関数として表されることを数式２４、数式２５は示している。ここで用いられるMonin-Obukov長Ｌｍｏは、以下の式で定義される。
【００３５】
【数２８】
【数２９】
ρは空気密度、Ｔｖは仮温度であり、Ｈ、λＥは地表面における顕熱および潜熱フラックスである。数式２６と数式２７のｓ＝θおよびｓ＝ｑに関する式をそれぞれｕ＊、Ｈ、λＥについて解くと、
【００３６】
【数３０】
【数３１】
【数３２】
で、数式３０〜数式３２は、無次元消散率の関数形が既知であれば、消散率を与えることで連立方程式として解くことができ、フラックスｕ＊、Ｈ、λＥが求められる。大気安定度パラメータζによる無次元関数表現は、接地境界層におけるMonin-Obukov similarity（ＭＯＳ）が成立することに基づいており、物理量の鉛直プロファイルがζで一意に表現できることを示している。例えば数式２４、数式２５中のφｍ≡ｋkarｚ／ｕ＊∂ｕ／∂ｚ、
φｓ≡ｋkarｚｕ＊／ｗ’ｓ’∂ｓ／∂ｚ をζの関数で表すのにBusinger-Deyer式がよく用いられる。また、無次元消散率の関数形も観測結果を元にいくつか提案されている。例えば、Kader(1992)による経験式は以下である。
【００３７】
【数３３】
【数３４】
温位と比湿の大気安定度に対する応答が相似であれば φεｃ＝φε０＝φεｑ が成立し、数式３０、数式３１、数式３２は繰り返し計算によって解くことが出来る。更にスカラーの無次元消散率の相似性の成立を仮定した場合は、ＣＯ２の消散率から、数式３１、数式３２と同様にしてフラックスを算出することが出来る。
【００３８】
次に、渦相関センサによる消散率の算定について説明する。シンチロメータ１との比較に用いるために、渦相関センサによる計測値から消散率を算定する方法について述べる。消散率の算定には、スペクトル（エネルギーや物質量変動の周波数表現）の慣性小領域における乱流特性を利用して間接的に求める慣性消散法（inertial dissipation method）を用いる。ただし、フーリエ変換の結果は一般に乱れが大きく、消散率の推定誤差要因となりうるため、ここではスペクトルの実空間表現である構造関数にもとづいた消散率の推定法を用いる。ＴＫＥおよびスカラー変動のスペクトルに対応する２次の構造関数Ｄab(r)は次式で定義される。
【数３５】
ここでａ，ｂはそれぞれｕ，θ，ｑなどの風速やスカラーを表し、ｒは実空間における距離である。時系列のデータ間の時間差をτとすると平均風速を用いて、τｕ＝ｒの関係にある。Kolmogorovの理論によればｒが慣性小領域にあるとき２次の構造関数と消散率は以下の関係にある。
【数３６】
【数３７】
ここでεとεｓはそれぞれＴＫＥおよびスカラー変動の消散率を表す。αuuはKolmogorov定数（およそ0.55）で、αssはObukhov-Corrsin定数（およそ0.8）と呼ばれる普遍的な定数である。比例定数Ｃs２は構造関数定数で、次のａとｂの組み合わせ［Ｔ：Ｔ］（気温）、［ｑ：ｑ］（比湿）および［Ｔ：ｑ］はそれぞれ構造関数定数ＣＴ２、Ｃｑ２およびＣＴｑに対応する。数式３６、数式３７に示すように２次の構造関数がｒ2/3に比例していることから、慣性小領域に相当するｒの範囲で構造関数にｒ2/3の曲線をフィットさせることで数式３６、数式３７から消散率を求め、消散法によるフラックス計算へと進む（図１１参照）。超音波風速計のパス長ｄｓより小さい乱流変動は測定されず、測定高さよりも大きい渦はその一部分しか測られないことから、十分に乱流変動計測が保障されるフィット範囲として ２ｄｓ＜ｒ＜（ｚ−ｄ０）／２ を選択する。
【００３９】
次に、シンチロメータ１を用いた消散率の算定方法について述べる。SLS40A（Scintec AG，Germany）は、市販のdisplaced-beam small aperture scintillometer で、波長0.67μｍのHeNeレーザーを偏光面が直交する２本のビームに分離し、２光路のシンチレーションからＴＫＥおよび気温変動の消散率を求め、フラックスを算出する装置である。SLSシリーズは単一光源を用いることで動作および計測の安定性と経済性を両立させた製品である。乱流による屈折率変動は大気中を伝搬する単色光の光強度のlog分散として計測することが出来る。独立した２光源と受光機がある場合の光強度のlog共分散Ｂ12は次式で表される（Thiemann，1992）。
【００４０】
【数３８】
ここで、ｘは、パス長Ｌの光軸上での位置；Ｋ＝２π／λは光の波数；ｋ乱流の波数；Ｃｎ２は屈折率変動の構造関数定数、ｌ０内部スケール（乱流渦が徐々に小さな渦に崩壊していく際、これ以上小さくなると分子粘性により渦が維持できなくなり熱として分解していくスケールの目安）である。Ｊ０およびＪ１は０次と１次の第１種ベッセル関数で、ビーム間の距離ｄおよび受光機の直径Ｄによって定まる感度分布を表す。検出器単体による分散を表すＢ１とＢ２は数式３８でｄ＝０とすることで得られる。数式３８中のφｎは屈折率変動の次元スペクトルで、
【数３９】
であり、ｆ（ｋｌ０）は、Hill and Clifford(1978)による屈折率変動の３次元スペクトルモデルである。数式３８は、σｘ２＜0.3の低散乱条件（Rytov approximation）が前提である。式中の被積分関数はシンチロメータ１の感度分布を表す荷重関数で、パスの中央で最大値となる釣鐘型の関数形を取る。数式３９を数式３８に代入し、ｘおよびｋについて積分することで次式が得られる。
【００４１】
【数４０】
ｆＢはｌ０，ｄおよびＤの増加に対するＢ12の減衰を表す関数である。相関係数ｒ12＝Ｂ12／（Ｂ１Ｂ２）１／２ ＝Ｂ12／Ｂ1＝Ｂ12／Ｂ2 はｌ０，ｄおよびＤの関数となり、ｄとＤを定数と考えると、相関係数ｒ12から内部スケールｌ０が求められる。ｌ０がわかれば、Ｂ1（あるいはＢ２）からＣｎ２が求められる。Scintecのソフトウェアでは、ｌ０に対するｒおよびｆＢをパス長Ｌに対する数表として用意し、Ｂ12とＢ1およびＢ2 からリアルタイムでｌ０とＣｎ２を計算する。得られたｌ０から乱流運動エネルギー（ＴＫＥ）の消散率は次の関係によって求められる。
【００４２】
【数４１】
ここでν（ｍ2ｓ-1）は空気の動粘性係数である。屈折率の変動に最も寄与するのは気温Ｔ（Ｋ）であり、ついで比湿ｑ（kgkg-1）の変動である。ＣＯ２などのガスは大気中に微量しか存在しないため屈折率変動には影響しない。屈折率変動の構造関数定数Ｃｎ２は気温変動の構造関数定数ＣＴ２（Ｋ２ｍ−２／３）、比湿変動Ｃｑ２（ｍ−２／３）および気温と比湿Ｔ−ｑの結合構造関数定数ＣＴｑ（Ｋｍ−２／３）によって次のように表される（Hill et al.，1980）。
【００４３】
【数４２】
定数ＡＴとＡｑは光の波長、気圧、気温および比湿の影響を受ける。本実施形態で使用するシンチロメータ１で利用する波長λ＝0.67μｍにおいては、Andreas(1989)による以下の式を用いる。
【数４３】
ここでＰ（ｈＰａ）は気圧である。光源に単一波長を用いるシンチロメータ１ではＣｎ２の値が一つしか得られないため数式４２からＣＴ２およびＣｑ２を一意に求めることが出来ない。この問題に対し、Thiermannは比湿がＣｎ２に与える影響が気温に比べて１オーダー小さいためこれを無視できるとし、数式４２の右辺第１項のみからＣＴ２を求めた。本実施形態で用いるシンチロメータ１に付属のフラックス計算ソフトウェアはこの計算方法を採用しており、運動量フラックスおよび顕熱フラックスのみが計算される。
【００４４】
また、本実施形態では水蒸気やＣＯ２密度などのスカラーについてもシンチロメータ１に対応する値を得るための手法を適用する。
【００４５】
まず水蒸気補正のための既存手法であるが、数式４２から明らかなように、シンチロメータ１は水蒸気の変動に対しても感度を持ち、水蒸気フラックスが顕熱フラックスに比べて非常に大きい条件ではＣＴ２が過大評価されてしまうことから水蒸気に対する補正の必要性が指摘されている（Green and Hayashi，1998）。水蒸気影響の古典的な補正法はWesely(1976)によって提唱されたボーエン比（顕熱と潜熱の比：Ｈ／λＥ）補正である。ボーエン比補正にはｗ’Ｔ’とｗ’ｑ’が必要なため、渦相関システムなどのセンサによる測定値を用いる。Moene(2003)はボーエン比補正について精査し、気温と比湿の標準偏差σＴ，σｑの比を用いる手法が、フラックスの比を用いた手法に対して、測定項目が少ないため誤差が小さいことを示している。すなわち
【数４４】
とすることで、数式４２は、Ｃｎ２とＣＴ２に関して次のように変形できる。
【数４５】
Ｔとｑの平均値と分散、そして相関係数ＲＴｑを他のセンサで計測することで、水蒸気の影響をキャンセルしたＣＴ２を計算することが出来る。
【００４６】
次にスカラーフラックスへの拡張であるが、水蒸気補正の際に導入した水蒸気と気温の構造関数定数比と分散比の相似性が成立するなら数式４４にしたがって、シンチロメータ１によるＣＴ２と、渦相関センサによるσＴ，σｑからＣｑ２を算出できると考えられる。同様にこの相似性をＣＯ２などのガスフラックスまで拡張すると以下が成り立つ。
【数４６】
これにより、渦相関センサによって得られた気温、水蒸気密度およびＣＯ２密度の分散とシンチロメータ１によって得られた気温変動の構造関数定数から水蒸気およびＣＯ２変動の構造関数定数ＣＴｑ２，Ｃｃ２が求められる。この算出法は純粋なシンチロメータ１による導出ではなく、渦相関センサとの組み合わせによることを留意する必要がある。しかし、現時点でシンチロメータ１と同様のスケールでの変動を測定する手法は存在しない。数式３７にｑおよびｃに関する構造関数定数を代入して消散率εｑとεｃが求められ、ε，εＴを加えて消散法によるフラックス計算へと進む。
【００４７】
続いて、渦相関法の補正法の誘導について説明する。
【００４８】
以下では、シンチロメータ１によるフラックスおよび消散率にはSASの添え字を、ECに相当する値にはECの添え字を添えて区別することとする。ここで、数式１６〜数式１８において各フラックスの比をとると数式４７〜数式４９となり、シンチロメータ１と渦相関センサによる消散率の比（εｉSAS／εｉEC）と渦相関法によるフラックスより、シンチロメータ１の空間スケールに相当するフラックスが得られる。なお、εｉにおける添字ｉは、何の消散率を表すかによってＴ，ｑ，ｃのいずれかの添字が該当することを意味している。
【数４７】
【数４８】
【数４９】
【００４９】
ここで、εＴSAS／εＴECは、乱流強度（水平風速で標準化した垂直風偏差：σｗ／ｕ）および風向に応答して変化することが分かっている。また、シンチロメータ１のパスに対する風向の変化はシンチロメータ１の空間スケールに影響する。つまり、風向の変化、すなわちパスに対する風の迎え角が変化することで、空間スケール（測定する範囲）が変化する。ただし、このほかに風速や大気の安定度なども影響するが、空間スケールの変化に最も大きく寄与するのが風向と考えられる。そこで、εｉSAS／εｉEC を乱流強度とシンチロメータ１の空間スケールを変数とする関数で表すことで、渦相関センサによる計測値からシンチロメータ１の空間スケールに対応するフラックスを算出することが可能となる。
【００５０】
シンチロメータ１の空間スケールに対応するフラックス算出に関しては以下のとおりである。まず、フラックス測定センサによって検出される乱流信号は風上の領域を起源とし、その広さ（ソースエリア）は風速や大気安定度、そしてセンサの測定パス長に対して変化することになる。シンチロメータ１と渦相関の空間スケールの違いを定量的に解析するために、Kormann and Meixner(2001) のフットプリントモデルを応用してシンチロメータ１のソースエリアを算出すると以下のとおりとなる。すなわち、フットプリントＦ（ｘ）は風上方向にｘ（ｍ）離れた位置のフラックスへの寄与を表した関数で、風速ｕ（ｍｓ−１）、摩擦速度ｕ＊（ｍｓ−１）測定高さｚ（ｍ）、大気安定度パラメータζを変数とする。さらに、De Bruin et al.(2002) の計算例に従い、シンチロメータ１のソースエリア（ＡSAS）を、Ｆ（ｘ）とセンサのパス長で表される横風方向ｙ（ｍ）の感度関数Ｄ（ｘ，ｙ）との積が描く等値線に囲まれた面積から求める（図６参照）。ここで、Ｄ（ｘ，ｙ）はシンチロメータ１が信号を検出する感度がパスの中央部で最も高いことを表す関数で、本実施形態では式形の簡単さから、Wang et al.(1978) が示した感度関数Ｄ（ｘ，ｙ）
［数５０］
Ｄ（ｘ，ｙ）＝（ｙ／Ｌｐ（１−ｙ／Ｌｐ））５／６
を標準化して用いる。ただし、Ｌｐはシンチロメータ１のパス長である。また、風向によって変化するシンチロメータ１のソースエリアは、パスに直交する方角からの風向角（θ）を用いて次のように表すこととする。
【数５１】
シンチロメータ１および渦相関センサのソースエリアの比ＲＡ（＝ＡSAS／ＡEC ）は、風向がシンチロメータ１のパスに平行な場合（θ＝９０°、２７０°）に１となり、直交風向時に増加する（図７参照）。
【００５１】
次に、領域平均フラックスを算出する。ここでは、まず上記で求めたソースエリア比ＲＡおよび乱流強度σｗ／ｕを変数とする関数ｆ（σｗ／ｕ，ＲＡ）の形を定め、一定期間に得られたεＴSAS／εＴEC の実測値に対する回帰によって、関数の係数を決定する。
【数５２】
次に、妥当と考えられるフラックス評価値を導くＲＡを求め、その関係が全ての風向に対して成立すると仮定して、数式４７より顕熱フラックスＨSASを求める。こうして求めたフラックスを領域平均フラックスと定義する。領域平均フラックスは全風向に対してシンチロメータ１の空間スケールを持ち、渦相関センサによる計測値のみでフラックス評価を行うことができるため、シンチロメータ１が測定不能な気象条件のフラックスも評価できる。さらに、水蒸気とＣＯ２が熱と同じ原理で輸送されると考えられることから、
【数５３】
と仮定すると、潜熱とＣＯ２についても同様に領域平均フラックスが計算できる（図８参照）。
【００５２】
なお、上述の実施形態は本発明の好適な実施の一例ではあるがこれに限定されるものではなく、本発明の要旨を逸脱しない範囲において種々変形実施可能である。
【実施例１】
【００５３】
中部地方の浅間東麓の５４年生の落葉広葉樹林においてＣＯ２を含むフラックス測定を、2002年から2004年を通して行った。ここでは、シンチロメータ１の空間スケールから決定した関数を用いて渦相関法の測定値を補正することにより、領域平均フラックスを算出するという「領域平均フラックスの算出」を実施し、算出した領域平均フラックスによって熱収支が正しく評価されていることを確認した。また、光合成量の計測結果とＣＯ２フラックスから求めた群落光合成量とを比較した。優先種であるダケカンバ個葉の光合成特性と葉の分布をもとに算出した森林の光合成ポテンシャルに対し、渦相関法が約ＣＯ２程度過小評価する一方で、領域平均フラックスによる算出値は良好に一致したことから、ＣＯ２交換量においても領域平均フラックスによって評価値の信頼性が向上することが分かった。さらに、森林のＣＯ２吸収量の評価を実施した結果、気象条件や樹種が類似した場所で得られた既往の渦相関法による評価結果とほぼ同等であった。以上につき、以下に実施例１として説明する。
【００５４】
ここでは、高さ２８ｍの鋼製足場材組上げ式のタワーを東西に８６ｍ離して配置し、東側の塔頂で渦相関法を実施するとともに領域平均フラックス測定のために、２基のタワー間をシンチロメータ１の測定パスとした。渦相関法では顕熱、潜熱、ＣＯ２フラックスを測定した。フラックス算出には１０Ｈｚサンプリング時系列の３０分ブロック平均を用いた。ＮＥＥの算出には降雨などの影響を受けにくく、測定値のドリフトが少ない閉光路（closed-path）式ＣＯ２分析計（LI-7000, Licor USA）によって計測したＣＯ２密度変動を用いた。内径４ｍｍのデカボンチューブ４０ｍを経由してタワー基部の観測舎に塔頂の大気サンプルを9.0Ｌmin−１で吸引し、その一部を2.0Ｌmin−１でＣＯ２分析計に導入した。サンプルエアが塔頂から分析器までに到達する時間遅れは、９秒程度であった。サンプルエアーの吸引によって信号の高周波成分が劣化してしまう問題は、乱流のスケールが大きな森林上ではほとんど考慮しなくて良いことが報告されている。当試験林においても、応答の速い開光路（open9-path）式のＣＯ２分析計と比較して、closed-path式のシステムが問題ない精度でＣＯ２フラックスを計測することを確認した。領域平均フラックスを計算するための関数ｆ（σｗ／ｕ，ＲＡ）を定めるために、シンチロメータ１の測定値が得られた2002年６月から１１月までの晴天日の測定値を用いて回帰計算を行った。
【００５５】
次に、欠測データの補間について説明する。フラックスからＮＥＥ（Net Ecosystem Exchange：生態系純交換量）を算出するためには、計測器の故障、メインテナンス、停電などさまざまな理由により発生する欠測を補い、連続した時系列を得る必要がある。そのために、1.5時間以内の短期欠測は線形補完を行い、それ以外はＣＯ２フラックスと相関の高い気象要素を用い、あらかじめ求めておいた気温‐呼吸量および光‐光合成の関係から補間を行った（図９参照）。測定値が得られない場合はAMeDASデータから当試験林の日平均日射量および気温を推定した。
【００５６】
また、フラックス測定値との比較を行うために、葉面積、個葉の光合成量の測定を行った。タワーを中心とした６０ｍ（Ｎ−Ｓ）×１４０ｍ（Ｅ−Ｗ）の範囲を１０ｍ×１０ｍの区画に分け、すべての生態データ取得の基準とした。ここでは、葉面積の季節変動を明らかにするために、開口部１ｍ×１ｍの落葉採取装置（リタートラップ）を林床から1.3ｍの高さに千鳥状に１９個配置し、９月中旬から１１月上旬の落葉期に落葉・落枝を定期的に採取した。また、強風などによる落葉期以外の落葉・落枝を考慮して着葉期の６月から９月上旬まで、平均的な落葉量を示した３個所のリタートラップを用いて落葉量を調査した。あらかじめ求めておいた落葉の乾燥重量と面積の１次回帰式から、採取した落葉サンプルの葉面積を求め、全葉面積から調査区域のＬＡＩ（Leaf Area Index：波面積指数）を求めた。光合成量の測定は、渦相関法を実施したタワーの北側に隣接したダケカンバ（Betula ermanii Cham.）の高度９ｍ付近にある任意の個葉を対象に、携帯型光合成蒸散量測定装置（LI-6400, Licor USA）を用いて行った。ダケカンバは、事前の毎木調査によって、当試験林を構成する樹木の地上部（枝＋幹）炭素量の44.3％を占める優先種であることが分かっている。当装置は小型のリーフチャンバを備え、チャンバ内の気温、湿度およびＣＯ２濃度を外部環境に追従させて制御できる。ダケカンバ葉が十分に光飽和する光強度から７段階（2000, 1500, 1000, 500, 200, 100, 0μmolm-2s-1）に変化させ、最も光合成能力が高いと思われる７月から９月期間で光強度に対する光合成量の反応（光−光合成反応）を調べた。
【００５７】
以上の結果、以下のような結果が得られた。まず、2002年から2004年にかけての年毎のフラックスデータ欠測率はそれぞれ１８％、５％、３％であった。2002年に欠測率が高いのは落雷による停電の影響で、そのほかにはメインテナンスや悪天候、機器の故障などの原因があった。ここで、FLUXNET(Falge et al.,2002)や、Asia Flux(Asia Flux運営委員会、2003) などでは信頼性の高いフラックス測定のために、乱流データのクオリティコントロール（ＱＣ）が提言されている。ＱＣでは乱流の成立性が乏しい状況での渦相関法は信頼性が低いとして、摩擦速度ｕ＊や時系列データの定常性に基準を設け、基準に満たないデータを棄却して補間の対象とする。このため、標準的な森林フラックスサイトの欠測率は３０％程度と言われている。しかし、本実施例では、渦相関法の信頼性が低いとされる乱流の乏しい状態あるいは対流の大きな状態における計測値も含めてシンチロメータ１のデータと比較し、その結果に基づいて領域平均フラックスを算出するため、ＱＣは適用せずすべてのデータを用いた。
【００５８】
続いて、補正関数の決定を行った。ここでは、シンチロメータ１と渦相関センサによる温度変動の消散率比εＴSAS／εＴEC は、機械的な乱流強度（σｗ／ｕ）が大きな中立時では、ＥＣセンサによる消散率がＳＡＳを上回り、（σｗ／ｕ）が小さな安定時および不安定時ではＳＡＳの消散率がＥＣより大きくなった（図１０（ａ）、図１０（ｂ）参照）。消散率とフラックスの絶対値は数式４７のように比例関係にあるため、この結果は、対流が卓越する条件や成層した大気条件では渦相関センサがフラックスを過小評価するというKanda et al.(2004)の結果と一致した。また、シンチロメータ１のソースエリアが渦相関センサに比べて大きくなる（ＲＡが大きくなる）につれてεＴSAS／εＴEC の分散の上端が大きくなる傾向が見られた（図１０（ｃ）、図１０（ｄ）参照）。センサの空間スケールが大きくなることで、より広い領域を起源とする信号の平均値を検出できるため、空間的スケールの大きな大気の乱流現象を捉えることができ、結果的に上記のσｗ／ｕに対するεＴSAS／εＴEC の関係が強調されると考えられた。
【００５９】
そこで、図１０（ａ）および図１０（ｂ）のｆ（σｗ／ｕ，ＲＡ）が大きい領域の漸近値をオフセットと考え、指数減衰関数（数式５４、数式５５）によって上記の関係を表現することとした（図１０の右欄参照）。ここでは、数式中の係数ａ１，ａ２，ａ３を昼間のおよび夜間について最小二乗法によって決定した。夜間と昼間で係数が異なるのは、ソースエリアの算出に用いたフットプリントモデルが夜間と昼間で異なる関数で表されるためである。
［数５４］
ｆ（σｗ／ｕ，ＲＡ）＝ｙ０＋ａ１exp［−（σｗ／ｕ−ｘ０）／ａ２］
［数５５］
ｘ０＝（ＲＡ−１）／ａ３
【００６０】
関数ｆ（σｗ／ｕ，ＲＡ）と消散率比εＴSAS／εＴEC の比較を図１２に示す。昼間（図１２（ａ））に比べて夜間（図１２（ｂ））の分散が大きいのは、渦相関センサによる温度変動の消散率εＴEC の算出精度が乱流状態の乏しい状況で低下したことが原因と考えられた。
【００６１】
次に、領域平均フラックスが得られると考えられる最適なＲＡを計算するために、年間の熱フラックスの収支とＲＡの関係を調べた（図１３）。昼間の不安定時および夜間安定時共に、ＲＡの増加にしたがってインバランス率Ｉ
［数５６］
Ｉ＝（Ｈ＋λＥ−Ｒｎ）／Ｒｎ
が０に近づき、熱収支が釣り合う効果が見られた。特に、夜間安定時はその効果が強かった。潜熱と顕熱の他に林冠内の大気や植物体に蓄えられる熱や地中伝熱（ 林冠内貯熱）が熱収支に関与するが、その値は大きくても１０数Ｗｍ−２であることと、ＲＡが２０を超える測定値がほとんど無く、外挿による誤差の拡大を防ぐ観点から、インバランス率＝−0.3を示すＲＡを求めた（夜間10.8、昼間23.1）。これによってｆ（σｗ／ｕ，ＲＡ）がσｗ／ｕの変数として表され、渦相関センサの測定値に基づいて領域平均フラックスが得られた。
【００６２】
ここで、領域平均フラックスの特徴について記す。領域平均フラックスは、σｗ／ｕが小さな領域では渦相関法によるフラックスより大きな値を示し、σｗ／ｕが大きな領域では渦相関法より小さな値を示す。その結果、夜間の放出および昼間の吸収量ピークが大きく見積もられ、日変動の振幅が拡大する。図１４に熱およびＣＯ２フラックスの測定例を示す。熱フラックスでは領域平均フラックスによる熱収支の改善が見られる一方で、大気が中立な条件では、領域平均フラックスが渦相関法を下回る状況も見られた（図１４のＡ参照）。また、熱収支の改善効果が小さい状況も見られた。林冠内貯熱の寄与が大きい時間帯が存在する可能性もあると考えられた。ＦＣについても同様に放出・吸収ピークが拡大する傾向が見られ（図１４のＢ参照）、ＦＣに基づいて計算される群落光合成量や、群落呼吸量も変化することになる。
【００６３】
さらに、領域平均フラックスと土壌呼吸量との比較を行った。土壌中の動物および根の呼吸と微生物による有機物の分解によってＣＯ２が放出される土壌呼吸と、地上部の植物体による呼吸が合算して森林からのＣＯ２放出となり、その放出量を群落呼吸量と呼ぶ。森林上で観測したＣＯ２からは、群落呼吸量に相当する情報が得られる。一方、土壌呼吸量は、チャンバ法などによって実測が可能である。そこで、ＣＯ２フラックスの妥当性を検証するために、群落呼吸量と土壌呼吸量を比較した。
【００６４】
まず、群落光合成量Ａｇ、群落呼吸量ＲｅとＮＥＥの間には
［数５７］
ＮＥＥ＝−Ａｇ＋Ｒｅ
の関係がある。夜間はＡｇ＝０なので、ＮＥＥが群落呼吸量となる。群落呼吸量は温度の相関が高く、温度に対する群落呼吸量のパラメタリゼーションには数式５８が用いられる（Lloyd and Taylor,1969）。
【数５８】
【００６５】
ここで、Ｔ(Ｋ)は温度、Ｔ０(Ｋ)は任意に決定する基準温度で、Ｒ０はＴ０(Ｋ)におけるＣＯ２フラックスを表す。また、Ｑ１０は温度係数で、基準温度におけるＣＯ２フラックスに対する、基準温度から１０Ｋ上昇したときのＣＯ２フラックスの比を表す。気温と群落呼吸量の関係（2004年）によると（図１５参照）、見易さのため σｗ／ｕ＞0.3以上のデータを示した。そこで、当該森林の群落呼吸量を解析するために、ＳＡＳおよびＥＣから求められる夜間のＣＯ２フラックス（ＦＣ）のデータと温度の関係から、パラメータＱ１０，Ｒ０を決定した（表１参照）。ここでは、冬季間に積雪することを考慮して、Ｔには地温ではなく地上高さ３ｍの位置の気温を用いた。なお、群落呼吸量Ｒｅは、土壌呼吸量（土壌中の生物による呼吸と根の呼吸）および植物地上部の呼吸量を含む。
【表１】
【００６６】
渦相関法にもとづく群落呼吸量のパラメータＱ１０と領域平均フラックスによるＱ１０は年間の差異において異なる傾向を示したが、いずれも平均気温との間には明確な傾向は見られなかった。次に、夜間の温度−呼吸量の関係が昼間でも同様に成り立つと仮定して、昼間の気温と数式５８から群落呼吸量を推定した。こうして求めた温度と群落呼吸量の関係と後述する群落光合成量とを用いて、数式５７から欠測の補間を行った。以前、土壌呼吸チャンバを用いて評価した当試験林の2002年における年間の土壌呼吸量は6.6tCha−１yr−１であったのに対し、渦相関法と領域平均フラックスにもとづく群落呼吸量はそれぞれ、7.4tCha−１yr−１と7.3tCha−１yr−１で（表２参照）、フラックス測定手法間に明確な相違は見られなかった。一方、2003年および2004年では、渦相関法の3.9tCha−１yr−１および7.2tCha−１yr−１に対して、領域平均フラックスは6.7tCha−１yr−１および8.9tCha−１yr−１と大きな値を示した。群落呼吸量は、土壌呼吸量に植物地上部の呼吸量が加わるため、フラックスによる評価値がチャンバー法を上回るのは妥当な結果といえる。ただし、温度に対する群落呼吸量の分散が大きく（図１５参照）、カーブフィッティングの際の誤差が大きいこととなった。領域平均フラックスから算出した群落呼吸量は土壌呼吸チャンバの測定結果と一致し、領域平均フラックスによって正しく群落呼吸量が評価されることが確かめられた。
【表２】
【００６７】
続いて、領域平均フラックスと群落光合成量との比較を行った。本実施例では、ＣＯ２フラックスの妥当性を検証するために、群落光合成量との比較を行った。ＣＯ２フラックスによって評価される群落光合成量は、森林全体の光合成量を表す。一方、葉の光合成量はリーフチャンバによって実測が可能である。しかし、森林は複雑な立体構造を持ち、林床面積あたりにその数倍の面積の葉が存在する。このため、個葉において計測された光合成量を群落光合成量と比較するためには個葉の計測結果を群落に拡張する必要がある。群落光合成量は、群落最上端に位置する葉の最大光合成量と葉面積から評価することができる。
【００６８】
まず第一に、群落光合成量のラメタリゼーションを行った。ＮＥＥと群落呼吸量の差は同化作用による群落光合成量（数式５９中のＡｇ）を示す。群落光合成量（μmolm-2s-1）Ｐｍは最大光合成量（μmolm-2s-1）を表し、αは初期光合成効率を表す。光合成の単位にはＣＯ２の分子量を乗じた密度表現（mgＣＯ２m-2s-1）も用いられる。数式５７に従って入射光量と数式５８に昼間の気温を代入して求めた呼吸量を差し引いたＦＣのプロットに数式５９をフィットしてパラメータを決定した。こうして求めたパラメータと入射光量および気温からＡｇとＲｅを算出し、欠測の補間を行った。
【数５９】
【００６９】
次に、個葉の光合成量を群落光合成量Ｐｃａｎに換算するための葉面積の算出結果を示す。落葉量から積算した葉面積指数（ＬＡＩ（Leaf Area Index：波面積指数））は、2002年から2004年までそれぞれ5.6，6.3，6.8で、増加傾向を示した。
【００７０】
また、図１６に示した葉面積指数の変化を見ると、展葉が完了した６月から、落葉が開始する９月までの期間も落葉が継続している。当試験林の優先種であるダケカンバは、春の展葉の後に２度目の伸張をすることが認められた。６月から９月の期間中も２次伸張による展葉が継続している場合、その間の落葉量をＬＡＩに算入するとＬＡＩを過大評価することになる。林床面積あたりに存在する実際の葉面積に近い値として、落葉開始直前の９月上旬の値を用いると、それぞれ、5.4，5.6，5.8となり、年間の格差は小さくなったが大きさの順位は変わらなかった。全落葉重量に占めるダケカンバの落葉重量は38.2％，42.1％，42.7％とほぼ41％程度であり、地上部炭素量の割合44.3％と近い値を示した。葉面積換算では29.1％，31.9％，32.5％とほぼ31％程度であった。
【００７１】
続いて、群落光合成量の算出を行った。領域平均フラックスによるＰｍは、７月および８月に最大値を示し、その値は1.7，1.9，1.8（mgＣＯ２m-2s-1）であった。一方、リーフチャンバによって計測したダケカンバ個葉の最大光合成量の平均値はそれぞれ0.8±0.2，0.9±0.1，0.9±0.1（mgＣＯ２m-2s-1）であった。リーフチャンバで得られたダケカンバ個葉の最大光合成量を群落最上端の最大光合成量とみなすと、Ｂｅｅｒの法則に従い、積算葉面積Ｆの位置で葉群に吸収された光量Ｉｉｎｋは以下で表される
【数６０】
ここで、Ｉ０は群落上端の入射光、Ｋは葉の角度などで決まる群落内の吸光定数で、0.5とする。ａは葉の透過率’（＝0.14）である。次に、入射光Ｉｉｎｋに対する光合成量Ｐは数式５９のPARをＩｉｎｋとすることで算出できる。図１７に、積算葉面積に対する葉群の吸光量Ｉｉｎｋおよび数式５９にＩｉｎｋを代入して求めた光合成量Ｐを示す。そして、ＰをＦで積分することにより、群落全体の光合成量Ｐｃａｎを計算することができる。
【数６１】
群落下層の植物は群落上端に比べて光合成能力が低い（最大光合成量が小さい）のが一般的で、実際にはさらに群落光合成量が小さくなることが考えられた。また、このモデルでは水分不足による気孔閉塞がもたらす昼寝現象などは考慮していない。したがって、本計算結果は群落光合成のポテンシャルと考えるのが適当である。植物体内の水分状態が気孔の開閉に影響したり、乱流の状態によってＣＯ２の交換量が変化するなど、フラックスの瞬時値はほぼ群落光合成のポテンシャル以下に分布すると考えられる。図１８に示すように、領域平均フラックスによるＡｇの分布の上端とＰｃａｎがほぼ重なる。このことは、領域平均フラックスが当該森林の光合成特性をよく表しているといえる。
【００７２】
また、図１９に当試験林で観測される最大光合成有効放射量2000μmolm-2s-1における群落光合成のポテンシャルとフラックスから導いた群落光合成量を示した。渦相関法の過小評価はＣＯ２フラックスにおいても明らかで、領域平均フラックスによる最大光合成量は、ＥＣによる最大光合成量の約1.5倍と、理論値である群落光合成のポテンシャルに近づいた。この結果、領域平均フラックスによって森林の光合成が正しく評価されることが示された。
【００７３】
続いて、領域平均フラックスの妥当性評価を行った。本発明者は、以前、シンチロメータ１の適用による熱収支式の整合をフラックス測定の妥当性の根拠とし、同じ原理で測定されるＣＯ２についても、フラックスが妥当に測定されると判断したが、本実施例においては更に領域平均フラックスから導いた群落呼吸量および群落光合成量と、群落呼吸量および個葉の光合成特性に基づく群落光合成量を比較した結果、領域平均フラックスによる算出値は、土壌呼吸量および光合成量の実測値とよく一致した。これにより、領域平均フラックスによるＣＯ２評価の妥当性が示された。また、熱フラックスの収支を判断材料として、フラックス測定の妥当性を評価する方法論が支持される結果となった。以上のような本実施例における結果は、当試験地特有の現象を捉えたものであるが、既存のフラックス測定サイトおよび新規測定サイトにおいても、渦相関法とシンチロメータ１と組み合わせた観測によってサイト固有の補正関数を求め、領域平均フラックスを求めることが可能と考えられた。
【００７４】
さらに、領域平均フラックスに基づく森林ＣＯ２吸収量の評価を行った。表２に示したように、ＮＥＥは大きな年変動を示した。これには、降水量、日射量、気温および水蒸気の飽差などの気象条件に加えて、窒素負荷や病害虫なども複合で影響すると考えられた。例えば、2004年は他の年より落葉時期が早く、９月中に葉量を減らしたにもかかわらず生産量は高かった。領域平均フラックスによるＮＥＥは渦相関法を上回り（図２０のＡ参照）、その結果、群落呼吸量および群落最大光合成量も大きく見積もられた（図２０のＢ参照）。３年間のＮＥＥの平均値は渦相関法では180gCm-2yr-1、領域平均フラックスでは266gCm-2yr-1と約1.5倍になった。以前、当試験林と類似した立地条件の４０年生のダケカンバ林で1998年から1999年にかけて実施された渦相関法によるフラックス測定では163gCm-2yr-1で、本実施例の結果に近い値であった。また、熱収支のインバランスを解消させるようにフラックスを補正したところ、214gCm-2yr-1となり約1.3倍となった。生態系呼吸量（Ｒｅ）と生態系純生産量（ＮＥＥ）の和で表されるＧＰＰ(Gross primary production)について比較すると、従来法で875gCm-2yr-1、補正後では1146gCm-2yr-1で、当試験林の３年間平均値811gCm-2yr-1および領域平均フラックスによる1156gCm-2yr-1と同等の結果であった。これは本実施例とは異なった補正法を用いているが、本実施例の結果と同様に渦相関法による過小評価が起こっていると考えられた。
【００７５】
以上のとおり、本実施例においては、森林におけるＣＯ２吸収量を高い精度で評価することを目的として、シンチロメータ１の空間スケールから決定した関数を用いて渦相関法の測定値を補正することにより、領域平均フラックスを算出した。ここでは、熱収支が整合するように調整した領域平均フラックスと、森林における群落呼吸量や群落光合成量との比較を通して領域平均フラックスの確かさを検証した。その結果、領域平均フラックスを基に算出した群落呼吸量は、以前に開発されたチャンバ法による土壌呼吸量の評価結果との一致を見た。また、従来法である渦相関法から算出した群落光合成量がダケカンバ個葉の計測結果から導いた群落光合成量を大きく下回るのに対し、領域平均フラックスによる群落光合成量は個葉の計測結果とよく対応することが明らかになった。以上より、森林における熱収支および群落のＣＯ２収支において、本発明者が領域平均フラックス評価法として開発した手法が従来の渦相関法の信頼性を高めることが分かった。例えば、従来法によって評価したＣＯ２吸収量の３年間の平均値は180gCm-2yr-1で、類似の立地条件にある落葉広葉樹林の値と同等であった。それに対し、本発明にかかる測定手法を用いて推定した森林によるＣＯ２吸収量の３年間の平均値は266gCm-2yr-1で、渦相関法による評価値の1.5倍であった。
【実施例２】
【００７６】
風速、乱流強度、ＴＫＥ（乱流運動エネルギー）、温度、比湿、ＣＯ２の各要素に関し、シンチロメータ１（ＳＣとも表記している）および渦相関センサ（ＳＣ）間の消散率比と風向の関係について実験を行った（図２１参照）。図中の角度は方位を表し、北を０°としている。グラフ中の値が大きければ大きいほどシンチロメータ１による消散率の比が大きく、例えば風速でいえば、北風ほどシンチロメータ１による消散率の比が大きくなる傾向が強い、という結果が得られた。
【実施例３】
【００７７】
摩擦速度ｕ＊、顕熱フラックスＨ、潜熱のフラックスλＥ、そしてＣＯ２フラックスｆｃのそれぞれに関し、シンチロメータ１および渦相関センサによるフラックス算出値の比較を行った（図２２参照）。図中、四角い記号で示すシンチロメータ１のパスに平行な西風、丸い記号でシンチロメータ１のパスに直交する北風のデータを示した。各期につき、シンチロメータ１による測定結果と渦相関センサによる測定結果の比は図に示すように、シンチロメータ１による測定結果が大きくなる傾向にあった。
【実施例４】
【００７８】
シンチロメータ１と渦相関センサのそれぞれの場合の非着葉期におけるインバランス率をグラフに表示することにより、熱収支インバランスが改善されたことの確認を行った（図２３参照）。インバランス率ＲＩ、熱収支式はそれぞれ図中に示す式で求めることができる。この結果、本発明にかかるＣＯ２フラックス測定法によれば従来法よりもインバランス率が改善されることがわかった（図２３参照）。
【図面の簡単な説明】
【００７９】
【図１】本発明のフラックス測定法におけるフラックスの定義と測定法を説明するための概略図である。
【図２】シンチロメータの測定原理を示す図である。
【図３】本発明におけるフラックス算出のデータ処理の流れを示す図である。
【図４】シンチロメータのパス長による測定領域の違いを渦相関センサの場合も併せて説明するための図である。
【図５】森林生態系のＣＯ２交換の概念を示す図である。
【図６】風向がパスに直交するときのフットプリントおよびソースエリアを示す図である。
【図７】ソースエリア比ＲＡと風向θの関係を示す図である。
【図８】領域平均フラックス算出の手順を示す図である。
【図９】欠測データの補間のフローを示す図である。
【図１０】εＴSAS／εＴEC に対する乱流強度σｗ／ｕおよびソースエリア比ＲＡの関係、ならびに補正関係ｆ（σｗ／ｕ，ＲＡ）と乱流強度σｗ／ｕおよびソースエリア比ＲＡの関係を示す図である。
【図１１】構造関数の例と曲線フィットによる消散率決定の内容（2002年６月10日）を示す図である。
【図１２】εＴSAS／εＴEC とｆ（σｗ／ｕ，ＲＡ）の比較を示す図で、ａ）は昼間、ｂ）は夜間を表している。
【図１３】ソースエリア比ＲＡとインバランス率Ｉ（ただし、Ｉは３年間の平均値）の関係を示す図である。
【図１４】熱フラックスおよびＣＯ２フラックスの日変化における渦相関法ＥＣと領域平均フラックスＳＡＳの関係を示す図である。
【図１５】気温と群落呼吸量の関係（2004年）を示す図であって、見やすさのためσｗ／ｕ＞0.3 の条件部分についてのみ表しているものである。
【図１６】リター積算に基づいた落葉期のＬＡＩの変化を示す図である。
【図１７】葉群の吸光量Ｉ０および光合成量（Ｐ）の積算葉面積（Ｆ）に対するプロファイルを示す図である。
【図１８】群落光合成のポテンシャルおよび群落光合成量の入射光量（PAR）との関係（2004年７月）を示す図である。
【図１９】群落光合成ポテンシャルと群落光合成量の比較（PAR＝2000μmolm-2s-1）を示す図である。
【図２０】領域平均フラックスと渦相関フラックスによる生態系純交換量（ＮＥＥ）および群落最大光合成量（Ａｇ）の比較を示す図である。
【図２１】シンチロメータ（ＳＣ）および渦相関センサ（ＥＣ）間の消散率比と風向の関係を各要素について示す図である。
【図２２】摩擦速度ｕ＊、顕熱フラックスＨ、潜熱のフラックスλＥ、ＣＯ２フラックスｆｃのそれぞれに関し、シンチロメータおよび渦相関センサによるフラックス算出値の比較結果を示す図である。
【図２３】シンチロメータと渦相関センサのそれぞれの場合の非着葉期におけるインバランス率が改善されたことを示す図である。
【符号の説明】
【００８０】
１ シンチロメータ
２ トランスミッタ
３ レシーバ
【特許請求の範囲】
【請求項１】
森林全体のＣＯ２吸収量を測定するためのＣＯ２フラックス測定法において、シンチロメータを用いた場合と渦相関法による場合の温度変動の消散率の比εＴSAS／εＴECが、シンチロメータを用いた場合の測定範囲ＡSASと渦相関センサによる場合の測定範囲ＡEC の比ＲＡ（＝ＡSAS／ＡEC ）との相関、および大気の乱流強度σｗ／ｕ との相関が高いことに基づき、これら２つの要素ＲＡ（＝ＡSAS／ＡEC ）およびσｗ／ｕ をパラメータとする関係式
【数１】
を、前記シンチロメータおよび前記渦相関センサを使って実際に得られた一定期間における測定結果に基づきあらかじめ用意しておくとともに、数式１中における測定範囲の比ＲＡには一定値を与えることによりすべての風向に対する前記シンチロメータの広い測定範囲を適用することとし、その後、森林の上空に実際に設置した渦相関センサによって大気の乱流強度σｗ／ｕを測定し、以下の関係式
【数２】
に基づき、渦相関法によって求めた顕熱フラックスＨECをシンチロメータによって求めた場合に相当する顕熱フラックスＨSASに補正し、尚かつ、渦相関法によって求めた潜熱フラックスλＥECおよびＣＯ２フラックスＦｃECに関しても、以下の関係式
【数３】
【数４】
に基づき、シンチロメータによって求めた場合に相当する潜熱フラックスλＥSASおよびＣＯ２フラックスＦｃSASに補正し、前記シンチロメータにて測定したのと同様の測定結果を得ることを特徴とする森林のＣＯ２フラックス測定法。
【請求項１】
森林全体のＣＯ２吸収量を測定するためのＣＯ２フラックス測定法において、シンチロメータを用いた場合と渦相関法による場合の温度変動の消散率の比εＴSAS／εＴECが、シンチロメータを用いた場合の測定範囲ＡSASと渦相関センサによる場合の測定範囲ＡEC の比ＲＡ（＝ＡSAS／ＡEC ）との相関、および大気の乱流強度σｗ／ｕ との相関が高いことに基づき、これら２つの要素ＲＡ（＝ＡSAS／ＡEC ）およびσｗ／ｕ をパラメータとする関係式
【数１】
を、前記シンチロメータおよび前記渦相関センサを使って実際に得られた一定期間における測定結果に基づきあらかじめ用意しておくとともに、数式１中における測定範囲の比ＲＡには一定値を与えることによりすべての風向に対する前記シンチロメータの広い測定範囲を適用することとし、その後、森林の上空に実際に設置した渦相関センサによって大気の乱流強度σｗ／ｕを測定し、以下の関係式
【数２】
に基づき、渦相関法によって求めた顕熱フラックスＨECをシンチロメータによって求めた場合に相当する顕熱フラックスＨSASに補正し、尚かつ、渦相関法によって求めた潜熱フラックスλＥECおよびＣＯ２フラックスＦｃECに関しても、以下の関係式
【数３】
【数４】
に基づき、シンチロメータによって求めた場合に相当する潜熱フラックスλＥSASおよびＣＯ２フラックスＦｃSASに補正し、前記シンチロメータにて測定したのと同様の測定結果を得ることを特徴とする森林のＣＯ２フラックス測定法。
【図６】
【図７】
【図１１】
【図１７】
【図２１】
【図２３】
【図１】
【図２】
【図３】
【図４】
【図５】
【図８】
【図９】
【図１０】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１８】
【図１９】
【図２０】
【図２２】
【図７】
【図１１】
【図１７】
【図２１】
【図２３】
【図１】
【図２】
【図３】
【図４】
【図５】
【図８】
【図９】
【図１０】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１８】
【図１９】
【図２０】
【図２２】
【公開番号】特開２０１１−１７７２８（Ｐ２０１１−１７７２８Ａ）
【公開日】平成２３年１月２７日（２０１１．１．２７）
【国際特許分類】
【出願番号】特願２０１０−２３２２８７（Ｐ２０１０−２３２２８７）
【出願日】平成２２年１０月１５日（２０１０．１０．１５）
【分割の表示】特願２００５−６３０９３（Ｐ２００５−６３０９３）の分割
【原出願日】平成１７年３月７日（２００５．３．７）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り ２００４年９月７日 日本農業気象学会、日本生物環境調節学会、日本植物工場学会、農業情報学会主催の「農業環境工学関連４学会２００４年合同大会」において文書をもって発表
【出願人】（０００１７３８０９）財団法人電力中央研究所 (1,040)
【Ｆターム（参考）】
【公開日】平成２３年１月２７日（２０１１．１．２７）
【国際特許分類】
【出願日】平成２２年１０月１５日（２０１０．１０．１５）
【分割の表示】特願２００５−６３０９３（Ｐ２００５−６３０９３）の分割
【原出願日】平成１７年３月７日（２００５．３．７）
【新規性喪失の例外の表示】特許法第３０条第１項適用申請有り ２００４年９月７日 日本農業気象学会、日本生物環境調節学会、日本植物工場学会、農業情報学会主催の「農業環境工学関連４学会２００４年合同大会」において文書をもって発表
【出願人】（０００１７３８０９）財団法人電力中央研究所 (1,040)
【Ｆターム（参考）】
[ Back to top ]