Climate Model Description

From MAGICC6 Wiki
Revision as of 11:16, 17 June 2013 by Antonius Golly (Talk | contribs) (Created page with "\subsection{From forcing to temperatures: the \\upwelling-diffusion climate model} \label{section_climatemodule} In the early stages, MAGICC's climate module evolved from the...")

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

\subsection{From forcing to temperatures: the \\upwelling-diffusion climate model} \label{section_climatemodule}

In the early stages, MAGICC's climate module evolved from the simple climate model introduced by \citet{Hoffert_1980_Role_DeapSea}. MAGICC's atmosphere has four boxes with zero heat capacity, one over land and one over ocean for each hemisphere. The atmospheric boxes over the ocean are coupled to the mixed layer of the ocean hemispheres, with a set of n-1 vertical layers below (see Fig.~\ref{fig_udebm_structure}). The heat exchange between the oceanic layers is driven by vertical diffusion and advection. In the previous model versions, the ocean area profile is uniform with depth and the corresponding downwelling is modeled as a stream of polar sinking water from the top mixed layer to the bottom layer. In this study, an updated upwelling-diffusion-entrainment (UDE) ocean model is implemented with a depth-dependent ocean area (from HadCM2). For simplicity, the following equations govern the uniform area upwelling-diffusion version of the model. 

Section~\ref{section_Appendix_implementation_UDE} provides details on the

UDE algorithms.




\subsubsection{Partitioning of feedbacks}

In order to improve the comparability between MAGICC and AOGCMs, and following earlier versions of MAGICC, we use different feedback parameters over land and ocean. This requires an adjustable land to ocean warming ratio in equilibrium based on AOGCM results. Given that in equilibrium the oceanic heat uptake is zero, the global energy balance equation can be written as:

\begin{eqnarray} \Delta Q_G=\lambda_G\Delta T_G=f_{L}\lambda_L \Delta T_{L} + f_{O}\lambda_O \Delta T_{O} \label{eq_globalenergybalance_equilibrium} \end{eqnarray} where $\Delta Q_G$, $\lambda_G$ and $\Delta T_G$ are the global-mean forcing, feedback, and temperature change, respectively. The right hand side uses the area fractions $f$, feedbacks $\lambda$, and mean temperature changes, $\Delta T$ for ocean ($O$) and land ($L$). As in earlier versions of MAGICC, the non-linear set of equations that determines $\lambda_O$ and $\lambda_L$ for a given set of equilibrium land-ocean warming ratio $RLO$ (=$\Delta T_L/\Delta T_O$), global-mean feedback $\lambda_G$, heat exchange and enhancement factors ($k$, $\mu$), is solved by an iterative procedure involving the set of linear Eqs.~(\ref{eq_fourboxequations_NO}--\ref{eq_fourboxequations_SL}), seeking the solution for $\lambda_L$ closest to $\lambda_G$. The procedure in version 6 has been modified slightly to take into account the time-constant radiative forcing pattern by CO$_2$ for the four boxes with hemispheric land/ocean regions, if prescribed.

%f9 \begin{figure}[t]\vspace*{2mm} \centering\includegraphics[width=8.3cm] {acpd-2007-0584-PartI-f09} 

\caption{The schematic structure of MAGICC's upwelling-diffusion energy balance module with
land and ocean boxes in each hemisphere. The processes for heat transport in the ocean are
deep-water formation, upwelling, diffusion, and heat exchange between the hemispheres. Not shown
is the entrainment and the vertically depth-dependent area of the ocean layers
(see Fig.~\ref{fig_udebm_tempprofile} and text).}

\label{fig_udebm_structure}\end{figure}



Following \citet{Wigley_Schlesinger_1985_AnalyticalSolutionsTemperature}, it is assumed that the atmosphere is in equilibrium with the underlying ocean mixed layer, so that the energy balance equation for the Northern Hemispheric ocean (NO) is: \begin{eqnarray} &f_{\rm NO}\lambda_{O}\Delta T_{\rm NO} = & \textrm{:infrared~outgoing~flux}\nonumber \\ &f_{\rm NO}\Delta Q_{\rm NO} & \textrm{:forcing} \nonumber \\ 

&+ k_{\rm LO}(\Delta T_{\rm NL} - \mu \Delta T_{\rm NO}) & \textrm{:land-ocean~heat~exchange} \nonumber \\
&+ k_{\rm NS}\alpha(\Delta T_{\rm SO} - \Delta T_{\rm NO}) & \textrm{:hemispheric

heat exch.}\label{eq_fourboxequations_NO} \end{eqnarray} where $\Delta T_{\rm NO}$ is the surface temperature change over the Northern Hemisphere ocean, $\Delta Q_{\rm NO}$ the radiative forcing over that region, $f_{\rm NO}$ the northern ocean's area fraction of the earth surface, $k_{\rm LO}$ the land-ocean heat exchange coefficient [W\,m{$^{-2}$}$^\circ$C$^{-1}$], a heat transport enhancement factor $\mu$ allowing for asymmetric heat exchange between land and ocean (1${\leq}\mu$ -- see Sect.~\ref{section_heatxchange_formulation} below), $k_{\rm NS}$ is the hemispheric heat exchange coefficient in the mixed layer. Following \citet{Raper_Cubasch_1996_Emulation_AOGCM_simplemodel} $\alpha$ is a sea-ice related adjustment factor to relate upper ocean temperature change to surface air temperature change (see Sect.~\ref{section_UD_equations}). Correspondingly, the equilibrium energy balance equations for the Northern Hemisphere land (NL), Southern Hemisphere ocean (SO) and Southern Hemisphere land (SL) are:

\begin{eqnarray} f_{\rm NL}\lambda_{L}\Delta T_{\rm NL} &=& f_{\rm NL}\Delta Q_{\rm NL} \nonumber \\ 

&&+ k_{\rm LO}(\mu \Delta T_{\rm NO} - \Delta T_{\rm NL})\label{eq_fourboxequations_NL}\\

f_{\rm SO}\lambda_{O}\Delta T_{\rm SO} &=& f_{\rm SO}\Delta Q_{\rm SO} \nonumber \\ &&+ k_{\rm LO}(\Delta T_{\rm SL} - \mu \Delta T_{\rm SO}) \nonumber \\ &&+ k_{\rm NS}\alpha(\Delta T_{\rm NO} - \Delta T_{\rm SO})\label{eq_fourboxequations_SO}\\ f_{\rm SL}\lambda_{L}\Delta T_{\rm SL} &=& f_{SL}\Delta Q_{SL} \nonumber \\ 

&&+ k_{\rm LO}(\mu \Delta T_{\rm SO} - \Delta T_{\rm SL})\label{eq_fourboxequations_SL}

\end{eqnarray}

As detailed below (Sect.~\ref{section_feedback_depending_forcings}), if the sensitivity factor $\xi$ is set different from zero (see Eq.~\ref{eq_feedback_dependency_forcing}), it is possible to make the feedback factors $\lambda$ in the energy balance equation dependent on the total radiative forcing. This forcing dependence of the feedback factors and the heat exchange enhancement factors are newly introduced in this version of MAGICC. The following two sections.~(\ref{section_heatxchange_formulation} and \ref{section_heatxchange_formulation}) are intended to provide both the motivation and details of these new parameterizations.


\subsubsection{Revised land-ocean heat exchange formulation} \label{section_heatxchange_formulation}


This section highlights a ``geometric effect that can cause effective climate sensitivities to change over time. The global-mean sensitivity may increase simply due to decreasing land-ocean warming ratios, given that climate feedbacks over land and ocean areas are different. To control the relative temperature changes over ocean and land, a heat transport enhancement factor $\mu$ is introduced. Enhancing the ocean-to-land heat transport ($\mu{\geq}$1) has the benefit that the simple climate model can better simulate some characteristic AOGCM responses. In the idealized forcing runs, AOGCMs often show a transient land-ocean warming ratio that slightly decreases over time, but stays above unity, combined with an increasing effective climate sensitivity in some models (see bottom rows in Fig.~\ref{fig_AOGCMcalibration_part1of3}, \ref{fig_AOGCMcalibration_part2of3}, 

and \ref{fig_AOGCMcalibration_part3of3}).

The higher land than ocean warming (RLO${>}$1) could be achieved by a smaller feedback (greater climate sensitivity) over land compared to the ocean boxes. However, as the land-ocean warming ratio decreases over time (due to less and less ocean heat uptake towards equilibrium), so would the effective global-mean climate sensitivity in previous model versions. The method used here, to allow both a RLO above unity and a non-decreasing effective climate sensitivity, assumes that ocean temperature perturbations influence the heat exchange more than land temperature changes. This asymmetric heat exchange formulation is then given by:

\begin{equation} {\rm HX}_{\rm LO}=k_{\rm LO}(\Delta T_{L} - \mu \Delta T_{O}) \label{eq_heatxchange} \end{equation} where HX$_{\rm LO}$ is the land-ocean heat exchange (positive in direction land to ocean), $\mu$ is the ocean-to-land enhancement factor and $\Delta T_L$ and $\Delta T_O$ are the temperature perturbations for the land and ocean region, respectively (cf.~Eq.~\ref{eq_fourboxequations_NO} ff.). Typical values for $\mu$ range between 1 and 1.4 as estimated from calibrating the CMIP3 ensemble (see Table~\ref{table_AOGCM_calibrationvaluesIII}).


\subsubsection{Accounting for climate-state dependent feedbacks} \label{section_feedback_depending_forcings}

Some AOGCM runs indicate higher effective climate sensitivities for higher forcings and/or temperatures. For example, the ECHAM5/MPI-OM model shows an effective climate sensitivity of approximately 3.5$^\circ$C after stabilization at twice pre-industrial CO$_2$ concentrations and 4$^\circ$C for stabilization at quadrupled pre-industrial CO$_2$ concentrations (see Fig.~\ref{fig_increasing_ClimSens_CCSM3_ECHAM5}b -- \citealp[see as well ][]{Raper_Gregory_Osborn_2001_diagnosingAOGCMresults, Hansen_etal_2005_Efficacies}). Given that the transient land-ocean warming ratio is the same for the 1pctto2$\times$ and 1pctto4$\times$ runs (see Fig.~\ref{fig_AOGCMcalibration_part1of3} last row), the 'geometric' effect discussed in the Sect.~\ref{section_heatxchange_formulation} would not explain this increase in climate sensitivity. An alternative explanation could be that climate feedbacks are climate-state dependent. The assumption in the standard energy balance Eq.~(\ref{eq_globalenergybalance}) with a constant global feedback ($\lambda$), with its attendant requirement that the outgoing energy flux scales proportionally with temperature change, may be an oversimplification. For example, the slow feedback due to retreating ice-sheets can lead to changes in the diagnosed effective sensitivities in AOGCMs (see e.g.~\citealp{Raper_Gregory_Osborn_2001_diagnosingAOGCMresults}) over long time-scales. \citet{Hansen_etal_2005_Efficacies} show that the 100-year climate response in the GISS model is more sensitive to higher forcings than to lower or negative forcings. Hansen et al.~(2005) express this effect by increasing efficacies for increasing radiative forcing. Table~1 in \citet{Hansen_etal_2005_Efficacies} suggests a gradient of roughly 1\% increase in efficacy for each additional Wm$^{-2}$ (OLS-regression of $E_a$ versus $F_a$ across the full range of CO$_2$ experiments), although some intervals (e.g.~from 1.25 to 1.5$\times$CO$_2$) show a slightly higher sensitivity of efficacy to forcing, i.e., 3\% per Wm$^{-2}$.

Rather than making the efficacies dependent on forcing, an alternative is to make the climate sensitivity dependent on the forcing level. This distinction, on whether to modify forcing or sensitivity, is not important when the climate system is at or close to equilibrium. However, if the efficacies of the forcing, instead of the feedback parameters are allowed to vary with forcing, the transient climate response after a change in forcing will be slightly faster. In this MAGICC version, if a forcing dependency of the sensitivity is assumed, the land and ocean feedback parameters $\lambda_L$ and $\lambda_O$ are scaled as

\begin{equation} \lambda = \frac{\Delta Q_{\rm 2\times}}{\frac{\Delta Q_{\rm 2\times}}{\lambda_{\rm 2\times}}+\xi(\Delta Q- \Delta Q_{\rm 2\times})} \label{eq_feedback_dependency_forcing} \end{equation} where $\lambda_{\rm 2\times}$ is the feedback parameter (=$\frac{\Delta Q_{\rm 2\times}}{\Delta T_{\rm 2\times}}$) at the forcing level for twice pre-industrial CO$_2$ concentrations. The sensitivity factor $\xi$ (KW$^{-1}$\,m$^{2}$) scales the climate sensitivity in proportion to the difference of forcing away from the model-specific ``twice pre-industrial CO$_2$ forcing level ($\Delta Q{-} \Delta Q_{\rm 2\times}$). The 1\% increase in efficacy for each additional unit forcing in Hansen's findings translates into a feedback sensitivity factor $\xi$ of 0.03\,KW$^{-1}$\,m$^{2}$ (assuming a climate sensitivity \emph{$\Delta T_{\rm 2\times}$} of 3\,$^\circ$C). Note that this scaling convention (Eq.~\ref{eq_feedback_dependency_forcing}) ensures that climate sensitivities are comparable for the equilibrium warming that corresponds to twice preindustrial CO$_2$ concentration levels (see Table~\ref{table_AOGCM_climsenscomparison}).


\subsubsection{Efficacies} \label{section_efficacies}

Efficacy is defined as the ratio of global-mean temperature response for a particular radiative forcing divided by the global-mean temperature response for the same amount of global-mean radiative forcing induced by CO$_2$ \citep[see Sect.~2.8.5 in][]{Forster_Ramaswamy_etal_2007_IPCCAR4_Chapter2_radiativeForcing}. In most cases, the efficacies are different for different forcing agents because of the geographical and vertical distributions of the forcing \citep{Boer_Yu_2003_ClimateSensitivityResponse, Joshi_etal_2003_improvedmetric_climatechange, Hansen_etal_2005_Efficacies}. The effective radiative forcing ($\Delta Q_e$) is the product of the standard climate forcing ($\Delta Q_a$), calculated after thermal adjustment of the stratosphere, and the efficacy (E$_a$). It is the effective forcings that are used in the energy balance equation (Eq.~\ref{eq_globalenergybalance}), although both effective and standard forcings are carried through in the MAGICC code. Note that this parameterization yields slightly faster transient climate responses compared to an approach where different climate sensitivities are applied for each individual forcing agent (cf.~Sect.~\ref{section_feedback_depending_forcings} above).

In MAGICC, forcings for some components differ by hemisphere and over land and ocean. Just as for the global sensitivity, this, in combination with different land/ocean feedback factors, results in MAGICC6 exhibiting efficacies different from unity for non-CO$_2$ forcing agents. In other words, efficacies different from unity are in part a consequence of the geometric effect described above. MAGICC calculates these internal efficacies using reference year (default 2005) forcing patterns. After normalizing these forcing patterns to a global-mean of \emph{$\Delta Q_{\rm 2\times}$} (default 3.71\,Wm$^{-2}$), the internal efficacy can be determined as

\begin{equation} E_{\rm int} = \frac{\Delta T_{\rm eff2\times}}{\Delta T_{\rm 2\times}}, \label{eq_efficacies_internal} \end{equation} where $T_{\rm eff2\times}$ is the actual global-mean equilibrium temperature change resulting from a normalized forcing pattern and \emph{$\Delta T_{\rm 2\times}$} is the corresponding warming for 2$\times$ CO$_2$ forcing, i.e., the climate sensitivity. For most forcing agents, these internal efficacies are very close to one, except for forcings with a strong land/ocean forcing contrast, such as aerosol forcings. For example, for direct aerosol forcing in the HadCM3 emulation (calibration III -- see Table~\ref{tableB3}) 

the efficacy is 1.14. By default, these internal efficacies are taken into

account when applying prescribed efficacies, so that:

\begin{equation} \Delta Q_e = \frac{E_a}{E_{\rm int}}\Delta Q_a \label{eq_efficacies_acc_for_intefficacies} \end{equation}

\newpage

\subsubsection{The upwelling-diffusion equations} \label{section_UD_equations}

The transient temperature change evolution is largely influenced by the climate system's inertia, which in turn depends on the nature of the heat uptake by the climate system. The transient energy balance equations can be written as:

\begin{eqnarray} f_{\rm NO}(\zeta_o \frac{d\Delta T_{\rm NO,1}}{dt}-\Delta Q_{\rm NO}+ \lambda_o \alpha \Delta T_{\rm NO,1} + F_{N}) =\nonumber\\ k_{\rm LO}(\Delta T_{\rm NL} - \mu\alpha\Delta T_{\rm NO,1})+k_{NS} \alpha(\Delta T_{\rm SO,1} - \Delta T_{\rm NO,1})\label{eq_transient_linear_equations_NO}\\ f_{\rm NL}(\zeta_L \frac{d\Delta T_{\rm NL}}{dt} - \Delta Q_{\rm NL} + \lambda_L \Delta T_{\rm NL}) =\nonumber\\ k_{\rm LO}(\mu\alpha\Delta T_{\rm NO,1}-\Delta T_{\rm NL})\label{eq_transient_linear_equations_NL}\\ f_{SO}(\zeta_o \frac{d\Delta T_{\rm SO,1}}{dt}-\Delta Q_{\rm SO}+ \lambda_o \alpha \Delta T_{\rm SO,1} + F_{S}) =\nonumber\\ k_{\rm LO}(\Delta T_{\rm SL} - \mu\alpha\Delta T_{\rm SO,1})+k_{\rm NS} \alpha(\Delta T_{\rm NO,1} - \Delta T_{\rm SO,1})\label{eq_transient_linear_equations_SO}\\ f_{\rm SL}(\zeta_L \frac{d\Delta T_{\rm SL}}{dt} - \Delta Q_{\rm SL} + \lambda_L \Delta T_{\rm SL}) =\nonumber\\ k_{\rm LO}(\mu\alpha\Delta T_{\rm SO,1}-\Delta T_{\rm SL})\label{eq_transient_linear_equations_SL} \end{eqnarray} where the adjustment factor $\alpha$ (default 1.2) determines -- over ocean areas -- the ratio of hemispheric changes in air ($\Delta T_{\rm xO}$) versus ocean mixed layer temperatures ($\Delta T_{\rm xO,1}$). Based on ECHAM1/LSG analysis \citep{Raper_Cubasch_1996_Emulation_AOGCM_simplemodel}, this sea-ice factor was first introduced by \citet{Raper_Gregory_Osborn_2001_diagnosingAOGCMresults} to account for the fact that the air temperature will exhibit additional warming, because the atmosphere feels warmer ocean surface temperatures where sea ice retreats. The bulk heat capacity of the mixed layer in each hemisphere $x$ is $f_x\zeta_o{=}f_x\rho c h_m$, where $\rho$ denotes the density of seawater (1.026$\times$10$^6$\,g\,m$^{-3}$), c is the specific heat capacity (0.9333\,cal\,g$^{-1}$$^\circ$C$^{-1}${=}4.1856$\times$0.9333\,Joule\,g$^{-1}$$^\circ$C$^{-1}$) and $h_m$ is the mixed layer's thickness [m]. The bulk heat capacity of the land areas is $f_x\zeta_L$, here assumed to be zero. The net heat flux into the ocean below the mixed layer is denoted by $F_x$.

Equation~(\ref{eq_transient_linear_equations_NL}) can then be written as: \begin{eqnarray} \Delta T_{\rm NL} = \frac{f_{\rm NL}\Delta Q_{\rm NL}+k_{\rm LO}\mu\alpha\Delta T_{\rm NO,1}}{f_{\rm NL}\lambda_L + k_{\rm LO}} \label{eq_transient_separating_TNL} \end{eqnarray} Substituting $\Delta T_{\rm NL}$ in Eq.~(\ref{eq_transient_linear_equations_NO}) yields: \begin{eqnarray} f_{\rm NO}(\zeta_o \frac{d\Delta T_{\rm NO,1}}{dt} - \Delta Q_{\rm NO} + \lambda_o \alpha\Delta T_{\rm NO,1} + F_N) = \nonumber \\ \frac{k_{\rm LO}}{\frac{k_{\rm LO}}{f_{\rm NL}}+\lambda_L}(\Delta Q_{\rm NL}-\lambda_L\mu\alpha \Delta T_{\rm NO,1})\nonumber\\ 

+ k_{\rm NS}\alpha(\Delta T_{\rm SO,1} - \Delta T_{\rm NO,1})

\label{eq_transient_TNL_into_TNO} \end{eqnarray} Provided we know the heat flux $F_N$ into the ocean below the mixed layer, we could now derive $d\Delta T_{\rm NO,1}/dt$. The net heat flux $F_N$ at the bottom of the mixed layer is determined by vertical heat diffusivity (diffusion coefficient $K_z$ [cm$^2$\,s$^{-1}${=}$3155.76^{-1}$\,m$^2$\,yr$^{-1}$]), and upwelling and downwelling (upwelling velocity $w$ [m yr$^{-1}$]), both acting on the perturbations $\Delta T$ from the initial temperature profile $T^0_{\rm NO,z}$. If the upwelling rate $w$ varies over time, the change in upwelling velocity $\Delta w^t{=}(w^t-w^0)$ compared to its initial state $w^0$ is assumed to act on the initial temperature profile, so that: \begin{eqnarray} F_N = \frac{K_z}{0.5h_d}\rho c (\Delta T_{\rm NO,1}-\Delta T_{\rm NO,2}) \nonumber \\ - w \rho c (\Delta T_{\rm NO,2} - \beta \Delta T_{\rm NO,1})\nonumber \\ - \Delta w \rho c (T^0_{\rm NO,2}- T^0_{\rm NO,sink}) \label{eq_transient_heatflux_bottom_mixedlayer} \end{eqnarray} where $T^0_{\rm NO,z}$ is the initial temperature for water in layer $z$ or in the downwelling pipe ($z$\,{=}\,``sink).

Given that the top layer is assumed to be mixed, the gradient of the temperature perturbations is calculated by the difference of the perturbations divided by half the thickness $h_d$ of the second layer (see Fig.~\ref{fig_udebm_tempprofile}). Substituting $F_N$ in Eq.~(\ref{eq_transient_TNL_into_TNO}) with Eq.~(\ref{eq_transient_heatflux_bottom_mixedlayer}) and transforming the equation to discrete time steps, yields: \begin{eqnarray} \frac{d\Delta T_{\rm NO,1}}{dt} \approx \frac{\Delta 

T_{\rm NO,1}^{t+1} - \Delta T_{\rm NO,1}^{t}}{\Delta t} =&  \\

\frac{1}{\zeta_o}\Delta Q_{\rm NO}^t & \textrm{:forcing} \nonumber \\ - \frac{\lambda_o \alpha}{\zeta_o}\Delta T_{\rm NO,1}^{t+1} & \textrm{:feedback}\nonumber \\ -\frac{K_z}{0.5h_d h_m}(\Delta T_{\rm NO,1}^{t+1}-\Delta T_{\rm NO,2}^{t+1}) & \textrm{:diffusion} \nonumber \\ + \frac{w^{t}}{h_m}(\Delta T^{t+1}_{\rm NO,2}-\beta \Delta T^{t+1}_{\rm NO,1}) & \textrm{:upwelling} \nonumber \\ + \frac{\Delta w^{t}}{h_m}(T^{0}_{\rm NO,2}-T^{0}_{\rm NO,sink}) & 

\textrm{:variable upwelling} \nonumber \\

+ \frac{k_{\rm LO}(\Delta Q^{t}_{\rm NL}-\lambda_L \mu\alpha\Delta 

T^{t+1}_{\rm NO,1})}{\zeta_o f_{\rm NO}(\frac{k_{\rm LO}}{f_{\rm
 NL}}+\lambda_L)} & \textrm{:land forcing} \nonumber \\

\nonumber + \frac{k_{\rm NS}\alpha}{\zeta_o f_{\rm NO}}(\Delta T_{\rm SO,1}^{t}-\Delta T_{\rm NO,1}^{t}) & \textrm{:inter-hemispheric ex.} \label{eq_transient_discrete_NO_mixed} \end{eqnarray} For the layers below the mixed layer (2$\leq$z$\leq$n--1), the temperature updating is governed 

by diffusion (first two

terms in Eq.~\ref{eq_transient_discrete_NO_submixed}) and upwelling (last two terms), so that: \begin{eqnarray} \frac{\Delta T^{t+1}_{{\rm NO},z}-\Delta T^{t}_{{\rm NO},z}}{\Delta t}= \nonumber \\ \frac{K_z}{0.5(h_d+h_d')h_d}(\Delta T^{t+1}_{{\rm NO},z-1}-\Delta T^{t+1}_{{\rm NO},z}) \nonumber \\ - \frac{K_z}{h_d^2}(\Delta T^{t+1}_{\rm NO,z}-\Delta T^{t+1}_{{\rm NO},z+1})\nonumber \\ + \frac{w^{t}}{h_d}(\Delta T^{t+1}_{{\rm NO},z+1}-\Delta T^{t+1}_{{\rm NO},z}) \nonumber \\ + \frac{\Delta w^{t}}{h_d}(T^0_{{\rm NO},z+1}-T^0_{{\rm NO},z}) \label{eq_transient_discrete_NO_submixed} \end{eqnarray} where $h_d'$ is zero for the layer below the mixed layer ($z${=}2) and $h_d$ otherwise, $\Delta w^t$ is the change from the initial upwelling rate.


For the bottom layer ($z = n$), the downwelling term has to be taken into account, so that: \begin{eqnarray} \frac{\Delta T^{t+1}_{{\rm NO},n}-\Delta T^{t}_{{\rm NO},n}}{\Delta t}= \frac{K_z}{h_d^2}(\Delta T^{t+1}_{{\rm NO,n}-1}-\Delta T^{t+1}_{{\rm NO},n}) \nonumber \\ + \frac{w^{t}}{h_d}(\beta\Delta T^{t}_{\rm NO,1}-\Delta T^{t+1}_{{\rm NO},n}) \nonumber \\ + \frac{\Delta w^{t}}{h_d}(T^0_{\rm NO,sink}-T^0_{{\rm NO},n}) \label{eq_transient_discrete_NO_bottom} \end{eqnarray}

Corresponding to the temperature calculations shown here for the Northern Hemisphere ocean (NO), the equivalent steps apply for the Southern Hemisphere ocean (SO). For simplicity, the equations described above are for the constant-depth area profile case, which MAGICC defaults to when the depth-dependency factor $\vartheta$ is set to zero. The detailed code for the general case with $0{\leq}\vartheta{\leq}1$ is given in Sect.~\ref{section_Appendix_implementation_UDE}.

\subsubsection{Calculating heat uptake} \label{section_Appendix_HeatUptake}

Heat uptake by the climate system can be calculated in different ways. One method is to use the global energy balance (Eq.~\ref{eq_globalenergybalance}). Using the effective sensitivity as in Eq.~(\ref{eq_globalenergybalance_equilibrium}) the heat uptake $F^t$ is estimated as: \begin{eqnarray} \frac{dH^t}{dt}=F^t = \Delta Q^t-(f_L \lambda_L \Delta T^t_{L} + f_O \lambda_O \Delta T^t_{O}) \label{eq_heatuptake_balance} \end{eqnarray}

%f10 \begin{figure}[t]\vspace*{2mm} \centering\includegraphics[width=8.3cm]{acpd-2007-0584-PartI-f10} 

\caption{The schematic oceanic area and initial temperature profiles in MAGICC's ocean hemispheres.
Diffusion driven heat transport is modeled proportional to the vertical gradient of temperature,
which is especially high below the mixed layer. } \label{fig_udebm_tempprofile}\end{figure}

For verification purposes MAGICC6 calculates heat uptake in two ways, both directly (as above) and by integrating heat content changes in each layer in the ocean (yielding identical results), given the assumed zero heat capacity of the atmosphere and land areas: \begin{eqnarray} \Delta H^t = \sum^{n}_{i=1} \frac{1}{\rho c h_i}\frac{(f_{\rm NO} \Delta T^t_{\rm NO,i} + f_{\rm SO} \Delta T^t_{\rm SO,i})}{f_O}+\epsilon \label{eq_heatuptake_integrating} \end{eqnarray} where $h_i$ is the thickness of the layer, i.e., $h_m$ for the mixed layer and $h_d$ for the others and $\epsilon$ is a small term to account for the heat content of the polar sinking water.


\subsubsection{Depth-dependent ocean with entrainment} \label{section_depthdependency_ocean}

\citet{Harvey_Schneider_1985_PartII, Harvey_Schneider_1985_PartI} introduced the upwelling-diffusion model with entrainment from the polar sinking water by varying the upwelling velocity $w$ with depth. Building on the work by \citet{Raper_Gregory_Osborn_2001_diagnosingAOGCMresults}, MAGICC6 also includes the option of a depth-dependent ocean area profile. If the depth-dependency parameter $\vartheta$ is set to 1 (default), a standard depth-dependent ocean area profile is assumed as in HadCM2 and used in \citet{Raper_Gregory_Osborn_2001_diagnosingAOGCMresults}. A constant upwelling velocity is assumed and mass conservation is maintained by ``entrainment from the downwelling pipe. With ocean area decreasing with depth and constant upwelling velocity, the upwelling mass flux would also have to decrease with depth. To offset this, the amount of entrainment into layer $z$ is assumed to be proportional to the decrease in area from the top to the bottom of each layer (cf.\ Fig.~\ref{fig_udebm_tempprofile}). We differ from the model structures tested by \citet{Raper_Gregory_Osborn_2001_diagnosingAOGCMresults}, by equating changes in the temperature of the entraining water to those in the downwelling pipe, namely a fraction $\beta$ (default 0.2) of the mixed layer temperature $\Delta T_{x,1}^{t-1}$ of the previous timestep in Hemisphere $x$. For a detailed description of the code, see the following Sect.~\ref{section_Appendix_implementation_UDE}. Simple upwelling-diffusion models can overestimate the ocean heat uptake for higher warming scenarios when applying parameter values calibrated to match heat uptake for lower warming scenarios (see e.g.~Fig.~17b in \citealp{Harvey_etal_1997_IPCC_IntroductionSimpleClimateModels}). To address this, MAGICC6 includes a warming-dependent vertical diffusivity gradient. The physical reasoning is that a strengthened thermal stratification and, hence, reduced vertical mixing leads to decreased heat uptake for higher warming. Thus, the effective vertical diffusivity at $K_{z,i}$ between ocean layer $i$ and $i$+1 is given by: \begin{eqnarray} K_{z,i} = {\rm max}\,(K_{z,{\rm min}}(1 - d_i)\frac{dK_{\rm z}}{dT}(\Delta T_{H,1}^{t-1}-\Delta T_{H,n}^{t-1})+K_z) \label{eq_verticalDiffusivity_heatdependent} 

\end{eqnarray}

where $K_{z,{\rm min}}$ is the minimum vertical diffusivity (default 0.1\,cm$^2$\,s$^{-1}$); $d_i$ is the relative depth of the layer boundary with zero at the bottom of the mixed layer and one for the top of the bottom layer; $\frac{dK_{\rm z}}{dT}$ is a newly introduced ocean stratification coefficient specifying how the vertical diffusivity $K_{\rm z}$ between the mixed layer 1 and layer 2 changes with a change in the temperature difference between the top/mixed and bottom ocean layer of the respective hemisphere at the previous timestep $t-$1 $(\Delta T_{H,1}^{t-1}{-}\Delta T_{H,n}^{t-1})$.


\subsection{Implementation of upwelling-diffusion-entrainment\\ equations} \label{section_Appendix_implementation_UDE}

This section details how the equations governing the upwelling-diffusion-entrainment (UDE) ocean (Eqs.~\ref{eq_transient_discrete_NO_mixed}, \ref{eq_transient_discrete_NO_submixed}, \ref{eq_transient_discrete_NO_bottom}) are implemented and modified by entrainment terms and depth-dependent ocean area (see Fig.~\ref{fig_udebm_tempprofile}). These equations represent the core of the UDE model and build on the initial work by \citet{Hoffert_1980_Role_DeapSea, Harvey_Schneider_1985_PartII, Harvey_Schneider_1985_PartI}.

The entrainment is here modeled so that the upwelling velocity in the main column is the same in each layer. Thus, the three area correction factors, $\theta_z^{\rm top}$, $\theta_z^{b}$ and $\theta_z^{\rm dif}$, applied below are:

\begin{eqnarray} \theta_z^{\rm top} = \frac{A_z}{(A_{z+1}+A_z)/2}\nonumber\\ \theta_z^{b} = \frac{A_{z+1}}{(A_{z+1}+A_z)/2} \nonumber\\ \theta_z^{\rm dif} = \frac{A_{z+1}-A_{z}}{(A_{z+1}+A_z)/2} \nonumber\\ \label{eq_areacorrection_thetatop} \end{eqnarray} where $A_z$ is the area at the top of layer $z$ or bottom of layer $z-1$ and the denominator is thus an approximation for the mean area of each ocean layer.

For the mixed layer, all terms in Eq.~(\ref{eq_transient_discrete_NO_mixed}) involving $\Delta T^{t+1}_{\rm NO,1}$ are collected on the left hand side in variable $A(1)$. All terms involving $\Delta T^{t+1}_{\rm NO,2}$ are collected in variable $B(1)$ on the left hand side. All other terms are held in variable $D(1)$ on the right hand side, so that the equation reads:

\begin{eqnarray} \Delta T_{\rm NO,1}^{t+1} = -\frac{B(1)}{A(1)}\Delta T_{\rm NO,2}^{t+1} + \frac{D(1)}{A(1)} \label{eq_udebm_coding_ALL1} \end{eqnarray} with \begin{equation}\label{eq_udebm_coding_A1} A(1) = 1.0+\theta_1^{\rm top}\Delta t\frac{ \lambda_O\alpha}{\zeta_o} \quad \textrm{:feedback over ocean} \end{equation} \[+\theta_1^{b}\Delta t\frac{ K_z}{0.5h_m h_d} \quad \textrm{:diffusion to layer 2} \] \[+\theta_1^{b}\Delta t\frac{ w^t \beta}{h_m} \quad \textrm{:downwelling} \] \[+\theta_1^{\rm top}\Delta t\frac{ k_{\rm LO}\lambda_L\mu\alpha }{\zeta_o f_{\rm NO} 

(\frac{k_{\rm LO}}{f_{\rm NL}} + \lambda_L)}  \quad  \textrm{:feedback over
land}\]

\begin{equation}\label{eq_udebm_coding_B1}B(1) = -\theta_1^{b}\Delta t\frac{ K_z}{0.5h_m h_d} \quad \textrm{:diffusion from layer 2} \end{equation} \[-\theta_1^{b}\Delta t\frac{ w^t}{h_m} \quad \textrm{:upwelling from layer 2}\] \begin{equation}\label{eq_udebm_coding_D1}D(1) = \Delta T_{\rm NO,1}^{t} \quad \textrm{:previous temp} \end{equation} \[+ \theta_1^{\rm top}\Delta t\frac{1}{\zeta_o}\Delta Q_{NO} \quad \textrm{:forcing ocean} \] \[+ \theta_1^{\rm top}\Delta t\frac{\alpha k_{NS}}{\zeta_o f_{NO}}(\Delta T^t_{\rm SO,1}-\Delta T^t_{NO,1}) \quad \textrm{:inter-hemis. exch.}\] \[+ \theta_1^{\rm top}\Delta t\frac{ k_{LO}\Delta Q_{NL}}{\zeta_o f_{NO} (\frac{k_{LO}}{f_{NL}} + \lambda_L)} \quad \textrm{:land forcing}\] \[+ \theta_1^{b}\Delta t\frac{\Delta w^t}{h_m}(T^0_{\rm NO,2}- T^0_{NO,sink}) \quad \textrm{:variable upwelling}\]


For the interior layers (2${\leq}z{\leq}n$), i.e., all layers except the top mixed layer and the bottom layer, the terms are re-ordered, so that A(z) comprises the terms for $\Delta T^{t+1}_{\rm NO,z-1}$, $B(z)$ the terms for $\Delta T^{t+1}_{\rm NO,z}$, C(z) the terms for $\Delta T^{t+1}_{\rm NO,z+1}$ and $D(z)$ the remaining terms, according to:

\begin{equation} \Delta T_{\rm NO,z-1}^{t+1} = -\frac{B(z)}{A(z)}\Delta T_{\rm NO,z}^{t+1} - \frac{C(z)}{A(z)}\Delta T_{\rm NO,z+1}^{t+1} + \frac{D(z)}{A(z)} \label{eq_udebm_coding_ALLmiddle} \end{equation} with \begin{equation} \scalebox{.84}[.84]{\lefteqn{A(z) = - \theta_z^{top}\Delta t\frac{K_z}{0.5(h_d+h_d')h_d} \quad \textrm{:diffusion from layer above} \label{eq_udebm_coding_Amiddle}}} \end{equation} \[B(z) = 1.0 + \theta_z^{b}\Delta t\frac{K_z}{h_d^2} \quad \textrm{:diffusion to layer below} \nonumber\] \[+\theta_z^{top}\Delta t\frac{K_z}{0.5(h_d+h_d')h_d} \quad \textrm{:diffusion to layer above} \] \begin{equation} +\theta_z^{top}\Delta t\frac{ w^t}{h_d} \quad \textrm{:upwelling to layer above} \label{eq_udebm_coding_Bmiddle} \end{equation} \[C(z) = - \theta_z^{b}\Delta t\frac{K_z}{h_d^2} \quad \textrm{:diffusion from layer below}\] \begin{equation} -\theta_z^{b}\Delta t\frac{ w^t}{h_d} \quad \textrm{:upwelling from layer below} \label{eq_udebm_coding_Cmiddle}\\ \end{equation} \[D(z) = \Delta T_{\rm NO,z}^{t} \quad \textrm{:previous temp}\] \[+\Delta t\frac{\Delta w^t}{h_d} (\theta_z^{b}T^{0}_{\rm NO,z+1}-\theta_z^{top}T^{0}_{\rm NO,z}) \quad \textrm{:variable upwelling}\] \[+\theta_z^{\rm dif}\Delta t\frac{ w^t}{h_d}\beta\Delta T_{\rm NO,1}^{t-1}\quad \textrm{:entrainment}\] \begin{equation} +\theta_z^{\rm dif}\Delta t\frac{\Delta w^t}{h_d} T^{0}_{\rm NO,sink} \quad \textrm{:variable entrainment} \label{eq_udebm_coding_Dmiddle} \end{equation} where $h_d'$ is zero for the layer below the mixed layer and $h_d$ otherwise. For the bottom layer, the respective sum factor $A(n)$ for $\Delta T^{t+1}_{\rm NO,n-1}$, $B(n)$ for $\Delta T^{t+1}_{\rm NO,n}$ and $D(n)$ for the remaining terms is:


\begin{eqnarray} \Delta T_{{\rm NO},n-1}^{t+1} = -\frac{B(n)}{A(n)}\Delta T_{{\rm NO},n}^{t+1} + \frac{D(n)}{A(n)} \label{eq_udebm_coding_Alln} \end{eqnarray} with \begin{equation} A(n) = - \theta_{n}^{\rm top}\Delta t\frac{K_z}{h_d^2} 

\quad  \textrm{:diffusion from layer n-1} \label{eq_udebm_coding_An}

\end{equation} \begin{equation}\label{eq_udebm_coding_Bn} B(n) = 1.0 + \theta_{n}^{\rm top}\Delta t\frac{K_z}{h_d^2} \quad \textrm{:diffusion to layer n-1}\end{equation} \[+\theta_{n}^{\rm top}\Delta t\frac{ w^t}{h_d} \quad 

\textrm{:upwelling to layer n-1} \]

\begin{equation}\label{eq_udebm_coding_Dn}D(n) = \Delta T_{{\rm NO},n}^{t} \quad \textrm{:previous temp} \end{equation} \[+\theta_{n}^{\rm top}\Delta t\frac{w^t}{h_d} \beta\Delta 

T^{t-1}_{\rm NO,1}  \quad  \textrm{:downwelling from top layer} \]

\[-\theta_{n}^{\rm top}\Delta t\frac{\Delta w^t}{h_d} T^{0}_{{\rm NO},n} \quad \textrm{:variable upwelling}\] \[+\theta_{n}^{\rm top}\Delta t\frac{\Delta w^t}{h_d} T^{0}_{\rm NO,sink} \quad \textrm{:variable downwelling} \]

With these Eqs.~(\ref{eq_udebm_coding_ALL1}--\ref{eq_udebm_coding_Dn}), the ocean temperatures can be solved consecutively from the bottom to the top layer at each time step.

Retrieved from "http://wiki.magicc.org/index.php?title=Climate_Model_Description&oldid=4"