飽和水蒸気圧の求め方 - もうカツ丼はいいよな

何年か前に飽和水蒸気圧の求め方に関して、Goff-Gratchの式やTetensの式やMurrayの式について少し触れたことがあった*1のだが、最近Wikipediaを見ていたらTetensの式だと思っていたものが

Tetens(1930)のパラメータ値によるAugust他の式
飽和水蒸気量 - Wikipedia

という表記になっていた。この部分はTetensの式となっていたものが割と最近更新された様子である。

Tetensの式はTetensの式ではない

Tetensの式やMurrayの式と呼ばれているこれらの式は、Magnusの式（後述するがこの呼称も適切ではない）と呼ばれている次式、
$e_s=C \exp \left( \frac{A t}{B + t} \right)$
のパラメータA, B, Cを定めたものに相当し、基本的には同様の式である。Tetensのパラメータがそうであるように、自然対数部分を常用対数として使うパラメータもある。
Tetens(1930)*2によるパラメータを用いたものがTetensの式、Murray(1966)*3のパラメータを用いたものがMurrayの式として呼ばれていたものに該当する。Tetensのパラメータが特に有名となっているのは、少ない桁数でそこそこ正確な近似を与えたという点で優れていたから、ということの様である。

August-Roche-Magnusの式

英語版 Wikipediaのクラウジウス-クラペイロンの関係式のページを見ると、飽和水蒸気圧の求め方が応用例としてあげられており、次のような記述がされている。

Fortunately, the August-Roche-Magnus formula provides a very good approximation, using pressure in hPa and temperature in Celsius:
$e_s(T)= 6.1094 \exp \left( \frac{17.625T}{T+243.04} \right)$ [7][8]
(This is also sometimes called the Magnus or Magnus-Tetens approximation, though this attribution is historically inaccurate.[9])
Clausius–Clapeyron relation - Wikipedia, the free encyclopedia

出典としてあげられているLawrence(2005)*4にもあるように、この式についての初期の仕事はAugust(1828)によってなされており、August-Rocheの式またはAugust-Roche-Magnusの式と呼ぶのがふさわしいようだ。このあたりに基づいて日本語版Wikipediaの記述も変更されたのだろう。

非August-Roche-Magnus式

August-Roche-Magnusの式の形をしていない飽和水蒸気圧を求める式として、Goff-Gratch(1946)の式や、Wexler(1976)の式、Sonntag(1990)の式などがあり、これらはより正確な飽和水蒸気圧を与えるとされている。
計算そのものはコンピュータを用いればすぐに終わるものだが、いずれも式が長く係数の桁数が多いので、そこそこの精度があれば良いという用途には不便である。

結局どうやって求めるの

英語版のクラウジウス-クラペイロンの関係式のページから引用した式の係数は、日本語の教科書によく載っている値や、日本語版Wikipediaに載っている値と若干違う。
この値は、Alduchov and Eskridge(1996)*5によりあたえられたもので、上に述べた3つの非Augst-Roche-Magnus式との相対誤差が、いずれに対しても0.4%以下となる（水上で-40℃〜50℃、氷上で-80℃〜0℃（氷上のものは係数が異なる）の範囲）ように選んだものである。有効な温度範囲も明示されているので、基本的にはこの式を使っておけば良いだろう。
Rの関数として定義する場合の例を示す。

AE1995 <- function(t, s.cool=FALSE){ # 入力：℃　出力：hPa
  ifelse(t > 0 | s.cool,
         6.1094 * exp(17.625 * t / (243.04 + t)), # 水上or過冷却水(s.cool=TRUE)
         6.1121 * exp(22.587 * t / (273.86 + t))  # 氷上
  )
}

Goff-Gratchの式との比較

Goff-Gratchの式をコントロールとして、AE1995を含むいくつかの式の相対誤差をプロットすると、以下のようになる。
f:id:Rion778:20160103135736p:plain
日常的な気温の範囲ではどれも大差なく、むしろAE1995は氷点以上のあてはまりがやや悪いようにも見えるが、AE1995はWexler(1976)やSonntag(1990)といった新しい式も含めた評価をしている点を注意しておく必要がある。面倒なので今回はやらないが、これらの式との比較をしてみるのも面白いだろう。
Okadaの式の当てはまりは特に良い様に見えるが、この式はGoff-Gratch式を最小二乗近似したもので、August-Roche-Magnusの式の一種ではない。そのため、元の式よりは簡単であるがやや入力が多い（スクリプトは後に示した）。
氷点以下の飽和水蒸気圧に関しては、AE1995がかなり良い値を示していることが分かるだろう。
なお、Tetensのパラメータは氷点以下のものがよく分からなかったので組み込んでおらず、氷点以下の値は適切ではない。後述のスクリプトを見てもらえばわかるが、Tetensのパラメータはとても簡単で、しかも常用対数を使っている。対象とする温度範囲に氷点以下が含まれないのであればこのパラメータを用いても問題ないだろう。この桁数でこの当てはまりは驚異的であり、Tetensの式と呼ばれてきただけの事はある。

Rのスクリプト

グラフの描画に用いたスクリプトを記しておく。
ライブラリは以下のものを使用した。

library(tidyr)
library(ggplot2)

飽和水蒸気圧を求める式を定義する。AE1995は前述の通り。

GofGra <- function(t, s.cool=FALSE){
  # Goff-Gratchの式
  ifelse((t > 0) | s.cool, 
    10^ # 水または過冷却水(s.cool=TRUE)の場合（-50〜100度）
      (10.79574 * (1 - 273.16/(t + 273.15)) -
         5.02800 * log((t + 273.15)/273.16, 10) +
         1.50475 * 10^{-4} * {1-10^(-8.2969 * ((t + 273.15)/273.16 - 1))} +
         0.42873 * 10^{-3} * {10^(4.76955*(1 - 273.16/(t + 273.15))) - 1} +
         0.78614), 
    10^ # 氷上の場合（-100度まで）
      (-9.09685 * (273.16/(t + 273.15) - 1) -
         3.56654 * (log(273.16/(t + 273.15), 10)) + 
         0.87682 * (1 - (t + 273.15)/273.16) + 
         0.78614)
    )
}

Okada <- function(t, s.cool=FALSE){
  # 岡田の式（Goff-Gratch式を最小二乗近似したもの）
  ifelse((t > 0) | s.cool,
    exp(1.809378 + # 水または過冷却水の場合（-30〜50度）
          0.07266115 * t +
          (-3.003879 * 10^-4) * t^2 +
          (1.181765 * 10^-6) * t^3 +
          (-3.863083 * 10^-9) * t^4),
    exp(1.809378 + # 氷上の場合（-30度まで）
          0.08238957 * t +
          (-2.990908 * 10^-4) * t^2 +
          (1.362765 * 10^-6) * t^3)
  )
}

# August式のパラメータ違い
Murray <- function(t){
  ifelse(t > 0,
         6.1078 * exp(17.2693882 * t /(t + 237.3)),
         6.1078 * exp(21.8745584 * t /(t + 265.5))
         )
  }
Tetens <- function(t){ 6.11 * 10^(7.5*t/(t + 237.3)) }

プロット。tidyrパッケージを使うとプロットに適した形式への変換が簡単だった。

# 対象とする温度範囲
ts <- seq(-100, 100, 0.1)
# 相対誤差を計算しデータフレームを作成
es.diff <- data.frame(air.temp = ts,
                      Okada =  (Okada(ts) - GofGra(ts)) / GofGra(ts) * 100,
                      Murray = (Murray(ts) - GofGra(ts)) / GofGra(ts) * 100,
                      Tetens = (Tetens(ts) - GofGra(ts)) / GofGra(ts) * 100,
                      AE1995 = (AE1995(ts) - GofGra(ts)) / GofGra(ts) * 100)
# longフォーマットへの変換(tidyr使用)
es.diff.l <- gather(es.diff, "method", "diff", -air.temp)
# plot
theme_set(theme_gray(base_family="HiraKakuProN-W3")) # 日本語プロット用設定
ggplot(es.diff.l, aes(x = air.temp, y = diff)) +
  geom_line(aes(color = method)) + xlim(-100, 100) + ylim(-1, 1) +
  labs(title = "Goff-Gratchの式の値に対する相対誤差",
       y = "相対誤差（％）", x = "気温（℃）")