%
%
%   2017/02/02  EBM
%               松田 幸樹
%
% 2017/12/22 作成
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%             Style  Setting             %%%%%%%%
% フォント: 12point (最大), 片面印刷

\documentclass[a4j,12pt,openbib,oneside,report,dvipdfmx]{jsbook}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%             Package Include            %%%%%%%%

\usepackage{Dennou6}            % 電脳スタイル ver 6
\usepackage{ascmac}
\usepackage{tabularx}
\usepackage{bmpsize}
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
%\usepackage{mathrsfs}
\usepackage{bm}
\usepackage{framed}
%\usepackage[dvipdfmx]{color}
\usepackage{docmute}
\usepackage{siunitx}
\usepackage{empheq}
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\usepackage{color}
\definecolor{shadecolor}{gray}{0.89}
\usepackage[hang,small,bf]{caption}
\usepackage[subrefformat=parens]{subcaption}
\captionsetup{compatibility=false}
\usepackage{wrapfig}
\usepackage{afterpage}
\usepackage{here}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%            PageStyle Setting           %%%%%%%%

\pagestyle{DAmyheadings}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%        Title and Auther Setting        %%%%%%%%
%%
%%  [ ] はヘッダに書き出される.
%%  { } は表題 (\maketitle) に書き出される.

\Dtitle{EBM 結果}   % 変更不可
\Dauthor{松田 幸樹}            % ゼミ担当者の名前
%\Ddate{2015/1/16}        % ゼミの日時 (毎回変更すること)
\Dfile{ebm\_conclusion.tex}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%   Set Counter (chapter, section etc. ) %%%%%%%%

\setcounter{chapter}{0}     % 章番号
\setcounter{section}{0}     % 節番号
\setcounter{subsection}{0}
\setcounter{equation}{0}   % 式番号
\setcounter{page}{0}      % 必ず開始ページは明記する
\setcounter{figure}{0}      % 図番号
\setcounter{table}{0}       % 表番号
%\setcounter{footnote}{0}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%        Counter Output Format           %%%%%%%%

\def\thechapter{\arabic{chapter}}
\def\thesection{\arabic{chapter}.\arabic{section}}
\def\thesubsection{\arabic{chapter}.\arabic{section}.\arabic{subsection}}
\def\theequation{\arabic{chapter}.\arabic{equation}}
\def\thepage{\arabic{page}}
\def\thefigure{\arabic{chapter}.\arabic{figure}}
\def\thetable{\arabic{chapter}.\arabic{section}.\arabic{table}}
\def\thefootnote{*\arabic{footnote}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%        Dennou-Style Definition         %%%%%%%%

%% 改段落時の空行設定
\Dparskip      % 改段落時に一行空行を入れる
%\Dnoparskip    % 改段落時に一行空行を入れない

%% 改段落時のインデント設定
\Dparindent    % 改段落時にインデントする
%\Dnoparindent  % 改段落時にインデントしない

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Macro defined by author

\def\univec#1{ \hat{ \Dvect{\rm #1}} }
%\def\DD#1#2{\frac{{\mathrm D} #1}{{\mathrm D} #2}}
%\def\dd#1#2{\frac{{\mathrm d} #1}{{\mathrm d} #2}}
\def\DD#1#2{\frac{D #1}{D #2}}
\def\dd#1#2{\frac{d #1}{d #2}}
\def\Dd#1{\; {\mathrm d} #1}
\def\D#1{\Dvect{#1}}
\def\cos#1{\mathrm{cos} #1}
\def\sin#1{\mathrm{sin} #1}
\def\tan#1{\mathrm{tan} #1}
\renewcommand{\DP}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
\def\tbref#1{表\ref{#1}}
\def\figref#1{図\ref{#1}}
\def\Del#1{\Delta #1}

%%%%%%%%             Text Start                 %%%%%%%%

\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{EBM を用いた太陽定数増減実験}

部分凍結解の初期値依存性を調べるために, 南北 1 次元エネルギーバランスモ
デル(EBM)を用いて実験を行う. 
EBM は, GCM よりも過程が単純であるため部分凍結解の初期値依存性の
解釈がしやすく, 計算量も少ないため様々な実験を行うことができる.
本章でははじめに, 作成した EBM について説明する. 
次に, 作成した EBM を用いて部分凍結解の初期値依存性について調べる.
ここでは, 水平格子点数や初期条件を変化させ, いくつかの実験を行う. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{EBM の概要}

本節では, 作成した EBM について詳細を説明する. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{支配方程式系}
\label{sec:ebm_dominant_eq}

本研究で作成したエネルギーバランスモデルは Sellers (1969) の定式化に従う.
モデルの支配方程式は, 地表面温度の時間発展の式
\begin{equation}
  C\DP{T\left(x,t\right)}{t} = F_{S} - F_{OLR} + \DP{}{x} D \left( 1 -
	 x^{2} \right) \DP{T\left(x,t\right)}{x},
	 \label{eq:time}
\end{equation}
で表される. ここで, $C$は熱容量, $T$は地表面温度, $x$は緯度の正弦, $D$は
拡散係数である. また, 右辺第一項$F_{S}$は太陽放射, 右辺第二項$F_{OLR}$は惑星放射,
右辺第三項は水平熱輸送を表す. ここで, 太陽放射$F_{S}$と惑星放射$F_{OLR}$
はそれぞれ
\begin{align}
 &F_{S} = Qs\left(x\right)\left(1-\alpha\left(x\right)\right) \label{solar} \\
 &F_{OLR} = A + BT(x,t) \label{olr}
\end{align}
である. ここで, $Q$は全球平均太陽放射フラックス, $\alpha$はアルベド,
$A,\ B$は定数である. 本研究では, Lindzen (1990) の値を参照し, $A=212.05\
\si{W.m^{-2}}$, $B=1.55\ \si{W.m^{-2}.K^{-1}}$とした. $s\left(x\right)$
は年平均太陽放射フラックスの子午面分布を表す関数であり,
\begin{equation}
 s\left(x\right) = 1 + s_{2}p_{2}\left(x\right) \label{s}
\end{equation}
と近似する. ここで, $s_{2}=-0.482$, $p_{2}\left(x\right)$は 2 次のルジャ
ンドル多項式で, $p_{2}\left(x\right)=\left(3 x^{2} - 1\right)/2$
である. $s_{2}$の値は Lindzen (1990) の値を参照した.  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{境界条件}

ここでは年平均の平衡状態について考えるため, 
計算領域を半球対称とし, 赤道($x=0$)から極($x=1$)までの領域を考える.
そして, 計算領域の両端において熱の出入りはないため, 境界条件は,
\begin{align}
 \left(\DP{T}{x}\right) &= 0,\ \ \ \left( x = 0 \right), \\
 \left(\DP{T}{x}\right) &= 0,\ \ \ \left( x = 1 \right)
\end{align}
である. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{支配方程式の離散化}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}[ht] 
 \centering
 \includegraphics[scale=0.5]{./Exp_EBM/fig/schematic_diagram.png}
 \caption{離散モデルにおける変数配置} \label{fig:mosiki}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ここでは, \ref{sec:ebm_dominant_eq}節に示した支配方程式の離散化について説明
する. \figref{fig:mosiki}に示すように, 赤道から極までを N 個の緯度帯に分割する. 
\eqref{eq:time}は下のように離散化する. 

\begin{align}
 C\frac{T^{l+1}_{k} - T^{l}_{k}}{\Delta t} =
  Q&s\left(x\right)\left(1-\alpha\left(x\right)\right)
  - \left(A + BT^{l}_{k}\right)
 \nonumber \\
 &+ \frac{D}{\Delta x^{2}}
  \left[ \left( 1 - x_{i+\frac{1}{2}}^{2}\right)\left(T_{i+1} - T_{i} \right)
  - \left( 1 - x_{i-\frac{1}{2}}^{2}\right)\left(T_{i} - T_{i-1} \right)
  \right].
\end{align}

ここで, 上付き添字$l$は時間を表し, 下付き添字$k$は緯度帯を表す. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{EBM における部分凍結解の初期値依存性} \label{sec:exp1}

本節では, GCM で見られるような部分凍結解の初期値依存性が EBM においても
見られるか否かを確認する. そのため, 解像度や初期条件を変化させ, いくつか
の実験を行う.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{実験設定}

アルベド$\alpha\left(x_{i}\right)$は,
\begin{align}
 \alpha\left(x_{i}\right) = \left\{
 \begin{array}{l}
  \alpha_{i} = 0.6\ \ \ \left(T\left(x\right) \le 271.15\ \si{K}\right), \\
  \alpha_{g} = 0.1\ \ \ \left(T\left(x\right) > 271.15\ \si{K}\right).
 \end{array}
 \right.
\end{align}
を与える. つまり, 凝結温度は$271.15\ \si{K}$とした.熱容量$C$は$1.0\
\si{J.K^{-1}}$とし, 拡散係数$D$は$0.2\ \si{W.m^{-2}.K^{-1}}$とする. モデ
ルの水平格子点数は$16$, $50$, $500$とする. ここで, 水平格子点数
$16$はGCM における T21 に対応する. 与えた全球平均太陽放射フラックスは,
$250\ \si{Wm^{-2}} \le Q \le 500\ \si{Wm^{-2}}$である.

各実験では次に述べる初期条件から平衡に達するまでモデルを積分することで平
衡状態を求める. \eqref{eq:time}の右辺の絶対値が$10^{-5}\ \si{W.m^{-2}}$
より小さくなると平衡状態になったと判断する.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{初期条件}

初期条件には次の 4 つの温度分布を与える(\tbref{tb:init_conf},
\figref{exp1:Init}). 一つ目は, 全緯度で$250\ \si{K}$である全球
凍結状態の温度分布(\figref{fig:Init_GN}, 青線), 二つ目は, $x< 0.40$にお
いて$300\ \si{K}$, $x \geq 0.40$において$250\ \si{K}$である部分凍結状態
の温度分布(\figref{fig:Init_PIC}, 緑線), 三つ目は, $x< 0.72$において$300\
\si{K}$, $x \geq 0.72$において$250\ \si{K}$である部分凍結状態の温度分布
(\figref{fig:Init_PIC}, 青線), 四つ目は, 全緯度で$300\ \si{K}$である氷な
し状態の温度分布である(\figref{fig:Init_GN}, 赤線). 以下では, 一つ目の初
期条件を用いた実験を [GIC250], 二つ目の初期条件を用いた実験を
[PIC6\_300-250], 三つ目の初期条件を用いた実験を [PIC11\_300-250], 四つ目の
初期条件を用いた実験を [NIC300] と呼ぶ(\tbref{tb:init_conf}).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 表

\begin{table}[ht]
 \begin{center}
  \caption{実験設定}
  \label{tb:init_conf}
  \begin{tabular}{|l|l|}
   \hline
   実験名 & 初期条件 \\
   \hline \hline
   [GIC250] & $T\left(x\right) = 250\ \si{K}\ \ (0\leq x \leq 1) $  \\
   \hline 
   [PIC6\_300-250] & $ T\left(x\right) =
       \begin{cases}
	300\ \si{K}\ \ (x< 0.40) \\
	250\ \si{K}\ \ (x \geq 0.40)
       \end{cases} $ \\
   \hline
   [PIC11\_300-250] & $ T\left(x\right) =
       \begin{cases}
	300\ \si{K}\ \ (x< 0.72) \\
	250\ \si{K}\ \ (x \geq 0.72)
       \end{cases} $ \\
   \hline
   [NIC300] & $T\left(x\right) = 300\ \si{K}\ \ (0\leq x \leq 1) $  \\
   \hline
  \end{tabular}
 \end{center}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 図

\begin{figure}[ht]
 
 \begin{tabular}{cc}
  \begin{minipage}{0.5\hsize}
   \centering
   \includegraphics[scale=0.4]{./Exp_EBM/fig/Init_GN.png}
   \subcaption{[GIC250] と [NIC300] に用いた初期温度分布. 赤線:
   [NIC300], 青線: [GIC250], 黒線: $271.15\ \si{K}$.}
   \label{fig:Init_GN}
  \end{minipage}
  \begin{minipage}{0.5\hsize}
   \centering
   \includegraphics[scale=0.38]{./Exp_EBM/fig/Init_PIC.png}
   \subcaption{[PIC6\_300-250] と [PIC11\_300-250] に用いた初期温
   度分布. 緑線: [PIC6\_300-250], 青線: [PIC11\_300-250], 黒線: $271.15\
   \si{K}$.}
   \label{fig:Init_PIC}
  \end{minipage}
 \end{tabular}
 \caption{水平格子点数 16 のときの初期温度分布}
 \label{exp1:Init}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{実験結果}

まず, 水平格子点数 16, $Q=300\ \si{W.m^{-2}}$として, [GIC250], [PIC6\_300-250],
[PIC11\_300-250], [NIC300] の実験を行った. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\figref{fig:time_exp1}に水平格子点数 16, $Q=300\ \si{W.m^{-2}}$のときの
地表面温度の時間発展を示す. \figref{fig:time_GIC}--\ref{fig:time_NIC}を
見ると, [GIC250], [PIC6\_300-250], [PIC11\_300-250], [NIC300] のいずれの
場合においても地表面温度はある解に収束している. [GIC250] は全球凍結状態,
[PIC6\_300-250], [PIC11\_300-250] は部分凍結状態, [NIC300] は氷なし状態
の平衡状態となっており, 多重解が得られた. さらに, [PIC6\_300-250],
[PIC11\_300-250] はどちらも部分凍結状態だが, 温度は異なっており, 結果と
して氷線も異なる緯度に位置する. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 図 時間発展

%全平均温度の方がよい？

 \begin{figure}[ht]

 \begin{tabular}{cc}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[scale=0.4]{./Exp_EBM/fig/Time_GIC.png}
   \subcaption{\hspace{0.25cm} [GIC250]}
   \label{fig:time_GIC}
  \end{minipage}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[scale=0.4]{./Exp_EBM/fig/Time_PIC6.png}
   \subcaption{\quad [PIC6\_300-250]}
   \label{fig:time_PIC5}
  \end{minipage}\\
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[scale=0.4]{./Exp_EBM/fig/Time_PIC11.png}
   \subcaption{\quad [PIC11\_300-250]}
   \label{fig:time_PIC11}
  \end{minipage}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[scale=0.4]{./Exp_EBM/fig/Time_NIC.png}
   \subcaption{\hspace{0.25cm} [NIC300]}
   \label{fig:time_NIC}
  \end{minipage}
 \end{tabular}
  \caption{水平格子点数 16, $Q=300\ \si{W.m^{-2}}$のときの地表面温度の時間発
  展(紫: $0\ \si{秒}$, 青: $5000\ \si{秒}$, 水色: $10000\ \si{秒}$,
  緑: $25000\ \si{秒}$, 黄: $50000\ \si{秒}$, オレンジ: $75000\ \si{秒}$,
  赤: 平衡時)}
  \label{fig:time_exp1}
 \end{figure}

% 色統一する
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
\figref{fig:Temp_5.11}に水平格子点数 16, $Q=300\ \si{W.m^{-2}}$のときの [PIC11\_300-250]
と [PIC6\_300-250] から得られた地表面温度の緯度分布を示す. 2 つの結
果を比較すると, 最も低緯度に位置する凍結した格子点は 3 格子点ずれている.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Temp-1

 \begin{figure}[ht]
  \centering
  \includegraphics[angle=-90,scale=0.33]{./Exp_EBM/fig/dcl_Temp-deg_16_6and11.pdf}
  \caption{水平格子点数 16, $Q=300\ \si{W.m^{-2}}$のときの平衡状態における地表面温度の緯度分
  布. 緑色の$\times$印は [PIC6\_300-250] の結果, 青色の$\times$印は
  [PIC11\_300-250] の結果, 黒線は凝固点$271,15\ \si{K}$である.}
  \label{fig:Temp_5.11}
 \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[PIC6\_300-250] と [PIC11\_300-250] において 最も低緯度に位置する凍結し
た格子点が 3 格子点ずれることを詳しく調べるために, [PIC5\_300-250] と
[PIC11\_300-250] だけでなく, 新たに 2 つの初期値を追加し
(\figref{fig:add_exp}の左図の赤線, オレンジ線)実験を行った.
\figref{fig:add_exp}の右図は\figref{fig:Temp_5.11}に追加実験の結果を追加
したものである. \figref{fig:add_exp}の右図より, 最も低緯度に位置する凍結
した格子点が 1 格子点ずつずれた部分凍結解が存在することが分かる.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Add Init and Temp-2

\begin{figure}[ht]
 
 \begin{minipage}{0.4\hsize}
  \centering
  \includegraphics[scale=0.38]{./Exp_EBM/fig/Init_PIC_Add.png}
 \end{minipage}
 \begin{minipage}{0.6\hsize}
  \centering
  \includegraphics[angle=-90,scale=0.3]{./Exp_EBM/fig/dcl_Temp-deg_16_6to11.pdf}
 \end{minipage}
 \caption{左図: [PIC6\_300-250] と [PIC11\_300-250] の初期温度分布と追加した初
 期温度分布(詳細は\tbref{init_add}参照). 
 右図: 水平格子点数 16, $Q=300\ \si{W.m^{-2}}$のときの平衡状態における地
 表面温度の緯度分布}
 \label{fig:add_exp}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 表

\begin{table}[ht]
 \begin{center}
  \caption{追加した初期温度分布の詳細}
  \label{init_add}
  \begin{tabular}{|c|l|}
   \hline
   色 & 初期条件 \\
   \hline \hline
   赤 & $ T\left(x\right) =
       \begin{cases}
	300\ \si{K}\ \ (x< 0.59) \\
	250\ \si{K}\ \ (x \geq 0.59)
       \end{cases} $ \\
   \hline
   オレンジ & $ T\left(x\right) =
       \begin{cases}
	300\ \si{K}\ \ (x< 0.66) \\
	250\ \si{K}\ \ (x \geq 0.66)
       \end{cases} $ \\
   \hline
  \end{tabular}
 \end{center}
\end{table}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

次に氷線緯度の太陽定数依存性を調べるために, $250\ \si{W.m^{-2}} \leq Q
\leq 500\ \si{W.m^{-2}}$の範囲で$Q$を$1\ \si{W.m^{-2}}$きざみで変化させ
て実験を行った. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\figref{fig:phase_diagram_16}に水平格子点数
16 のときの氷線緯度と太陽定数の関係を示す.
氷線緯度は, 凝固点の前後の格子点の温度の線形補間で求めた. 
全球凍結解, 部分凍結解, 氷なし解の 3 つの解や多重解が得られた.
また, 初期条件によって, 2 種類の部分凍結解が存在することが確認でき, 
[PIC11\_300-250] より [PIC6\_300-250] から得られた平衡解の氷線の方が低緯
度に位置している.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% フェーズダイアグラム

 \begin{figure}[ht]
  \centering
  \includegraphics[angle=-90,scale=0.35]{./Exp_EBM/fig/dcl_latitude_16.pdf}
  \caption{水平格子点数 16 のときの太陽定数と氷線緯度の関係. ピンクは
  [GIC250] の結果, 青は [PIC6\_300-250] の結果, 緑は [PIC11\_300-250] の
  結果, 赤は [NIC300] を示す. }
  \label{fig:phase_diagram_16}
 \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

最後に, 氷線緯度と太陽定数の関係の水平格子点数依存性について調べるため,
水平格子点数を 50, 100, 500 として同様の実験を行った.
\figref{fig:phase_exp1}は水平格子点数 16, 50, 100, 500 における氷線緯度
と太陽定数の関係である. \figref{fig:phase_exp1}より, 部分凍結解の初期値
依存性の幅は水平格子点数が増えると小さくなることがわかる. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% フェーズダイアグラム解像度依存性

 \begin{figure}[ht]

 \begin{tabular}{cc}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[angle=-90,scale=0.25]{./Exp_EBM/fig/dcl_latitude_16.pdf}
   \subcaption{水平格子点数 16}
   \label{fig:phase16}
  \end{minipage}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[angle=-90,scale=0.25]{./Exp_EBM/fig/dcl_latitude_50.pdf}
   \subcaption{水平格子点数 50}
   \label{fig:phase50}
  \end{minipage}\\
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[angle=-90,scale=0.25]{./Exp_EBM/fig/dcl_latitude_500.pdf}
   \subcaption{水平格子点数 500}
   \label{fig:phase500}
  \end{minipage}
 \end{tabular}
  \caption{太陽定数と氷線緯度の関係: (a) 水平格子点数 16, (b) 50, (c)
  500. ピンクは [GIC250] の結果, 青は [PIC6\_300-250] の結果,
  緑は [PIC11\_300-250] の結果, 赤は [NIC300] の結果を示す. (b), (c)
  には比較のために (a) の結果を黒点で描画している. }
  \label{fig:phase_exp1}
 \end{figure}

% 初期値統一しないとね

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\afterpage{\clearpage}
\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{まとめと考察}

本節では, GCM で見られた部分凍結解の初期値依存性が EBM でも見られるか調
べた. その結果, EBM においても部分凍結解の初期値依存性が存在する
ことや, その部分凍結解の初期値依存性の幅が水平格子点数を増やすと小さくなること
が分かった.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 各項の緯度分布

 \begin{figure}[H]
  \centering
  \includegraphics[scale=0.5]{./Exp_EBM/fig/Flux_PIC5.png}
  \caption{水平格子点数 16, $Q=300\ \si{W.m^{-2}}$のとき得られた
  [PIC6\_300-250] における平衡状態
  における太陽放射フラックス, 長波放射フラックス, 水平熱輸送
  フラックスの単位面積当たりの加熱率. 赤線は太陽放射フラックス, 青線は
  長波放射フラックス, 緑線は水平熱拡散フラックスの加熱率である}
  \label{exp1:flux}
 \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

従って, 水平格子点数を無限に大きくすると, 部分凍結解は一意に求まるように
思われる. つまり, 離散系における表現に部分凍結解の初期値依存性の原因があると想定される.
そこで, \eqref{eq:time}の右辺のうちどの項が特に離散化の影響を受けるのか
を調べることにする.
\figref{exp1:flux}に, 水平格子点数 16, $Q=300\ \si{W.m^{-2}}$のとき得ら
れた [PIC6\_300-250] 平衡状態における太陽放射フラックス, 長波放射フラックス, 水平熱輸送の
加熱率を示す. 
拡散型の水平熱輸送を用いたことで温度は連続的な緯度分布をしているため, 温度の一次関数である
長波放射フラックスは連続的な分布となっている. ところが, 吸収される太陽放
射フラックスと正味の水平熱輸送の加熱率は氷線近傍で不連続的な分布を示して
いる. 水平熱輸送の加熱率は他の二項と釣り合うように決定されるため, ここで
は太陽放射フラックスの加熱率に着目する. \eqref{solar}
と\eqref{s}より, $Qs(x)$は連続的な分布であるが, $\alpha$は不
連続な分布になっている. つまり, 氷線近傍における太陽放射フラックスの加熱
率が不連続な分布になるのはアルベドの分布に起因していると考えられる. そこ
で, 次節では離散モデルでのアルベドの表現を見直し, それによる解への影響を
調べることにする. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\afterpage{\clearpage}
\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{EBM におけるアルベドの表現の影響} \label{sec:exp2}

\ref{sec:exp1}節より, 部分凍結解の初期値依存性は支配方程式の離散化に原因があるこ
とがわかった. また, 太陽放射フラックスに着目すると氷線近傍において不連
続となっており, これはアルベドの分布に起因していると考えた. そこで本節で
は,  \ref{sec:exp1}節のアルベドの表現 (\figref{fig:albedo}の黒線) からサ
ブグリッドスケールの氷の面積を考慮し連続系における表現に近づけたアルベド
の表現 (\figref{fig:albedo}の赤線) に変更し, 追加実験を行う. そして, そ
の結果と\ref{sec:exp1}節の実験結果を比較する. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{実験設定}

本節では, サブグリッドスケールの氷の面積を考慮し, アルベド分布を考慮する.
このとき, アルベド$\alpha\left(x_{i}\right)$は, 
\begin{align}
 \alpha\left(x_{i}\right) &=
 \frac{1}{\Del{x}}\int_{\Del{x}}\alpha\left(x\right)dx \nonumber \\
 &= \frac{1}{\Del{x}}\int_{\Del{x}}\left[\alpha_{i}\beta\left(T\right) +
 \alpha_{g}\left(1 - \beta\left(T\right)\right)\right]dx
\end{align}
とする. ここで,
\begin{align}
 \beta\left(T\right) = \left\{
 \begin{array}{l}
  1\ \ \ \left(T\left(x\right) \le 271.15\ \si{K}\right), \\
  0\ \ \ \left(T\left(x\right) > 271.15\ \si{K}\right).
 \end{array}
 \right.
\end{align}
\begin{equation}
 T\left(x\right) = T\left(x_{i}\right) + \frac{T\left(x_{i+1}\right) -
  T\left(x_{i}\right)}{x_{i+1}-x_{i}}\left(x - x_{i}\right) \label{eq:temp}
\end{equation}
である. つまり, \eqref{eq:temp}よりアルベドは格子点の間の温度分布を線形と仮定し
て, 各格子内の温度が凝固点を下回る面積を推定し
(\figref{fig:albedo_menseki}の青い領域), 面積の重みをつけた平均によって
格子点のアルベドを求める. その結果,
例えば, 水平格子点数$16$, $Q=300\ \si{W.m^{-2}}$のときの平衡状態の温度分布
(\figref{fig:albedo_menseki}の赤線)におけるアルベド分布は, \figref{fig:albedo}の赤線の
ようになる. その他の設定は\ref{sec:exp1}と同じである. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% アルベド面積とアルベド分布

\begin{figure}[ht]
   \centering
   \includegraphics[scale=0.75]{./Exp_EBM/fig/albedo_menseki.png}
   \caption{水平格子点数$16$, $Q=300\ \si{W.m^{-2}}$のときの平衡状態にお
   ける地表面温度分布. 青で塗られた領域はある格子内の温度が凝固点を下回る領域で
   あり, 青線で囲まれ白で塗られた領域はある格子内の温度が凝固点を上回る領域である.}
   \label{fig:albedo_menseki}
\end{figure}

\begin{figure}[H]
 \centering
 \includegraphics[scale=0.75]{./Exp_EBM/fig/albedo.png}
 \caption{水平格子点数$16$, $Q=300\ \si{W.m^{-2}}$のときの平衡状態にお
 けるのアルベド分布. 黒線はアルベドの扱いを変更する前, 赤線はアルベド
 の扱いを変更した後のアルベド分布である. }
 \label{fig:albedo}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{初期条件}

実験名および初期条件は\ref{sec:exp1}と同じである (\tbref{tb:init_conf},
\figref{exp1:Init}). 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{実験結果}

まず, 水平格子点数 16, $Q=300\ \si{W.m^{-2}}$として, [GIC250], [PIC6\_300-250],
[PIC11\_300-250], [NIC300] の実験を行った. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\figref{fig:time_exp2}に水平格子点数 16, $Q=300\ \si{W.m^{-2}}$のときの
地表面温度の時間発展を示す. \figref{fig:time_GIC-2}--\ref{fig:time_NIC-2}を
見ると, [GIC250], [PIC6\_300-250], [PIC11\_300-250], [NIC300] のいずれの
場合においてもある解に収束している. [GIC250] は全球凍結状態,
[PIC6\_300-250], [PIC11\_300-250] は部分凍結状態, [NIC300] は氷なし状態
の平衡状態となっており, 多重解が得られた. さらに, [PIC6\_300-250],
[PIC11\_300-250] はどちらも部分凍結状態であり同じ温度となっており, 氷線
も同じ緯度に位置する. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 図 時間発展

 \begin{figure}[ht]

 \begin{tabular}{cc}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[scale=0.4]{./Exp_EBM/fig/Time_GIC2.png}
   \subcaption{\hspace{0.25cm} [GIC250]}
   \label{fig:time_GIC-2}
  \end{minipage}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[scale=0.38]{./Exp_EBM/fig/Time_PIC6-2.png}
   \subcaption{\quad [PIC6\_300-250]}
   \label{fig:time_PIC5-2}
  \end{minipage}\\
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[scale=0.4]{./Exp_EBM/fig/Time_PIC11-2.png}
   \subcaption{\quad [PIC11\_300-250]}
   \label{fig:time_PIC11-2}
  \end{minipage}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[scale=0.38]{./Exp_EBM/fig/Time_NIC2.png}
   \subcaption{\hspace{0.25cm} [NIC300]}
   \label{fig:time_NIC-2}
  \end{minipage}
 \end{tabular}
  \caption{水平格子点数 16, $Q=300\ \si{W.m^{-2}}$のときの地表面温度の時間発
  展(黒: $0\ \si{秒}$, 紫: $5000\ \si{秒}$, 青: $10000\ \si{秒}$, 水
  色: $50000\ \si{秒}$, 緑: $100000\ \si{秒}$, 黄: $150000\ \si{秒}$, オレ
  ンジ: $200000\ \si{秒}$, 赤: 平衡時)}
  \label{fig:time_exp2}
 \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

次に氷線緯度の太陽定数依存性を調べるために, $250\ \si{W.m^{-2}} \leq Q
\leq 500\ \si{W.m^{-2}}$の範囲で$Q$を$1\ \si{W.m^{-2}}$きざみで変化させ
て実験を行った. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\figref{fig:phase_diagram_16-2}に水平格子点数
16 のときの氷線緯度と太陽定数の関係を示す. 全球凍結解, 部分凍結解, 氷な
し解の 3 つの解や多重解が得られた. また, アルベドの表現を変更したことに
よって, 初期条件によらず部分凍結解は一意に求められている. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% フェーズダイアグラム

 \begin{figure}[ht]
  \centering
  \includegraphics[angle=-90,scale=0.4]{./Exp_EBM/fig/dcl_latitude_16-2.pdf}
  \caption{水平格子点数 16 のときの太陽定数と氷線緯度の関係. ピンクは
  [GIC250] の結果, 青は [PIC6\_300-250] の結果, 緑は [PIC11\_300-250] の
  結果, 赤は [NIC300] の結果を示す. }
  \label{fig:phase_diagram_16-2}
 \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

次に, 氷線緯度と太陽定数の関係の水平格子点数依存性について調べるために,
水平格子点数を 50, 100, 500 として同様の実験を行った.
\figref{fig:phase_exp2}は水平格子点数 16, 50, 100, 500 における氷線緯度
と太陽定数の関係である. これによると, 部分凍結解は水平格子点数によらずおおよそ
一意に求まることが分かる.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% フェーズダイアグラム解像度依存性

 \begin{figure}[ht]

 \begin{tabular}{cc}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[angle=-90,scale=0.25]{./Exp_EBM/fig/dcl_latitude_16-2.pdf}
   \subcaption{格子点数 16}
   \label{fig:phase16-2}
  \end{minipage}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[angle=-90,scale=0.25]{./Exp_EBM/fig/dcl_latitude_50-2.pdf}
   \subcaption{格子点数 50}
   \label{fig:phase50-2}
  \end{minipage}\\
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[angle=-90,scale=0.25]{./Exp_EBM/fig/dcl_latitude_500-2.pdf}
   \subcaption{格子点数 500}
   \label{fig:phase500-2}
  \end{minipage}
 \end{tabular}
  \caption{太陽定数と氷線緯度の関係: (a) 水平格子点数 16, (b) 50, (c)
  500. ピンクは [GIC250] の結果, 青は [PIC6\_300-250] の結果,
  緑は [PIC11\_300-250] の結果, 赤は [NIC300] の結果を示す. (b), (c) には比較のために (a) の結果を黒点で描画している. }
  \label{fig:phase_exp2}
 \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\figref{fig:hikaku}は各解像度ごとの氷線緯度と太陽定数の関係を
\ref{sec:exp1}節のアルベドの扱いを変更する前の実験と\ref{sec:exp2}節のアルベド
の扱いを変更した後の実験を比較したものである. \figref{fig:hikaku}より, アルベ
ドの扱い変更後の実験において得られる解はアルベドの扱い変更前の実験におい
て水平格子点数を大きくしたときに収束する解とおおよそ一致する. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 比較

 \begin{figure}[ht]

 \begin{tabular}{cc}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[angle=-90,scale=0.25]{./Exp_EBM/fig/dcl_latitude_16_1-2.pdf}
   \subcaption{格子点数 16}
   \label{fig:hikaku16}
  \end{minipage}
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[angle=-90,scale=0.25]{./Exp_EBM/fig/dcl_latitude_50_1-2.pdf}
   \subcaption{格子点数 50}
   \label{fig:hikaku50}
  \end{minipage}\\
  \begin{minipage}[b]{0.5\linewidth}
   \centering
   \includegraphics[angle=-90,scale=0.25]{./Exp_EBM/fig/dcl_latitude_500_1-2.pdf}
   \subcaption{格子点数 500}
   \label{fig:hikaku500}
  \end{minipage}
 \end{tabular}
  \caption{氷線緯度と太陽定数の関係の各解像度ごとの比較: (a) 水平格子点
  数 16, (b) 50, (c) 500. 赤はアルベドの扱いを変更する前, 青はアルベ
  ドの扱いを変更した後の結果である. }
  \label{fig:hikaku}
 \end{figure}

% 初期値統一しないとね

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\afterpage{\clearpage}
\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{まとめと考察}

本節では, サブグリッドスケールの氷の面積を考慮してアルベドの表現を変更し,
部分凍結解の初期値依存性について調べた. その結果, 部分凍結解の初期値依存性は
現れなかった. これは水平格子点数によらない結果である. また, アルベドの表現を変
更する前の実験において水平格子点数を大きくしたときに収束した結果とアルベ
ドの表現を変更した実験の結果はおおよそ一致した. 以上より, EBM における部分凍結解の
初期値依存性はサブグリッドスケールのアルベド分布を無視したことに起因して
いることがわかった. そして, 本研究で用いたアルベドの表現は上手く機能して
おり, それを用いることで部分凍結解の初期値依存性はなくなることがわかった. 

\ref{sec:exp1}節では, アルベドの表現を変更する前の実験において, 水平格子
点数を大きくすると部分凍結解の初期値依存性の幅は小さくなることがわかった.
従って, \ref{sec:exp1}節と同様のアルベドの表現をしている本研究の GCM に
おいても解像度を大きくすると, 部分凍結解の初期値依存性の幅は小さくなるこ
とが期待されるが, 容易ではない. しかし, \ref{sec:exp2}と同様に, GCM にお
いてもサブグリッドスケールの氷の面積を考慮したアルベド分布を与えることで第
\ref{exp_gcm}章で見られた部分凍結解の初期値依存性は小さくなると考えられ
る. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{初期値依存性の解釈}

ここでは, アルベドの表現を変更する前になぜ EBM において部分凍結解の初期値依
存性が現れたのかについて改めて考察することにする. それは, モデルの振る舞
いの理解に有益であるためである. 本節では,
\eqref{eq:time}の左辺をゼロとした平衡状態におけるエネルギーバランスの式,
\begin{equation}
 F_{S} - F_{OLR} + \DP{}{x} D \left( 1 - x^{2} \right)
  \DP{T\left(x,t\right)}{x}= 0 \label{eq:insta}
\end{equation}
を用いて, この式の解の振る舞いを考える. 任意の全球平均太陽放射$Q$に対し
て, 以下のように\eqref{eq:insta}を解く:
\begin{enumerate}
 \item 始めに氷線を仮定する. この時仮定した氷線を$x_{0}$とする.
 \item 氷線を仮定したことによりアルベド分布が求まり, \eqref{eq:insta}を
       解くことができる.
 \item \eqref{eq:insta}を解いて得られた温度分布から氷線を求める. この時
       求めた氷線を$x$とする.
 \item このとき, $ x = x_{0} $ならば$x$は\eqref{eq:insta}の解となる. 
\end{enumerate}
ここで, 物理定数の値やアルベドの与え方は\ref{sec:exp1}節と同様である. ただし, 上記の方
法では, 小氷冠不安定や大氷冠不安定のような不安定解も解として求まることに注
意しておく.

$Q=300\ \si{W.m^{-2}}$のときの仮定した$x_{0}$と$x$の関係を
\figref{fig:x0-x}に示す. 赤線は水平格子点数 16 のときの結果, 青線は水平
格子点数 500 のときの結果である. また, 緑線は$x = x_{0}$の線であり, 緑線
と赤線もしくは青線との交点が解となる. まず, 水平格子点数 16 のとき,
$x_{0}$に対応する$x$が階段状に分布する.
これは, 仮定する氷線の位置$x_{0}$を変化させても, $x_{0}$が隣の格子点の位
置に至るまではアルベドの分布に変化がなくそれに対応する$x$の値も変化しな
いためである. このため, 緑線と赤線の交点となっている部分凍結解が多数の
$x$で実現する. 一方, 水平格子点数 500 では一定の$x$となる$x_{0}$の範囲が
小さいので,
部分凍結解はおおよそ一意に求まることがわかる. 本節の実験では水平格子点数
500 を用いたが, 水平格子点数を無限にすれば連続系の表現と考えることができ,
このとき$x$は連続な分布となると推測される. 従って, 連続系では部分凍結解
は一意に求まると考えられる. また, \ref{sec:exp2}節のようにアルベドの表現
を変更すると, $x_{0}$の変化に対して連続的にアルベドの分布が変化するこ
とになる. これによって, \ref{sec:exp2}節では部分凍結解が一意に求められたと
考えられる.

GCM では, EBM には含まれていない多数の複雑な過程が取り入れられている. し
かし, 本節で述べた離散化された EBM の振る舞いと同様のことは GCM において
も起こっているだろう. GCM で見られた異なる氷線を持つ部分凍結解の存在には,
ここで述べた機構が重要な役割を果たしていると考えられる. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% x_0-x

 \begin{figure}[H]
  \centering
  \includegraphics[scale=0.6]{./Exp_EBM/fig/Interpretation.png}
  \caption{$Q=300\ \si{W.m^{-2}}$のときの$x_{0}$と$x$の関係 (青点は水平
  格子点数 16 のときのアルベドの扱い変更前の結果, 緑点は水平格子点数 500
  のときのアルベドの扱い変更前の結果, 緑点は水平格子点数 16 のときのアル
  ベドの扱い変更後の結果, 黒線は$x = x_{0}$の線).}
  \label{fig:x0-x}
 \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}

%%%%%%%%              Text End                  %%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%