\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small Tenth MSU Conference on Differential Equations and Computational Simulations. \newline \emph{Electronic Journal of Differential Equations}, Conference 23 (2016), pp. 213--221.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2016 Texas State University.} \vspace{9mm}} \begin{document} \setcounter{page}{213} \title[\hfilneg EJDE-2016/Conf/23 \hfil An improved method] {An improved method for estimating ice line for zonal energy balance climate models} \author[Z. Zhang, et al. \hfil EJDE-2016/conf/23 \hfilneg] {Zhenbu Zhang, Tor A. kwembe, Remata S. Reddy, Jerryl Roberts, Brittany Keys, Tsion Andine} \address{Zhenbu Zhang (corresponding author)\newline Department of Mathematics and Statistical Sciences\\ Jackson State University\\ Jackson, MS 39217, USA} \email{zhenbu.zhang@jsums.edu} \address{Tor A. kwembe \newline Department of Mathematics and Statistical Sciences\\ Jackson State University\\ Jackson, MS 39217, USA} \email{tor.a.kwembe@jsums.edu} \address{Remata S. Reddy \newline Department of Physics, Atmospheric Sciences \& Geoscience\\ Jackson State University\\ Jackson, MS 39217, USA} \email{remata.s.reddy@jsums.edu} \address{Jerryl Roberts \newline Department of Physics, Atmospheric Sciences \& Geoscience\\ Jackson State University\\ Jackson, MS 39217, USA} \email{jerryl.m.roberts@students.jsums.edu} \address{Brittany Keys \newline Department of Mathematics and Statistical Sciences\\ Jackson State University\\ Jackson, MS 39217, USA} \email{Brittanykeys.edu@gmail.com} \address{Tsion Andine \newline Department of Biology\\ Jackson State University\\ Jackson, MS 39217, USA} \email{tsion1995@yahoo.com} \thanks{Published March 21, 2016} \subjclass[2010]{35B30, 35K20, 35K57} \keywords{Temperature profile; climate change; energy balance; ice line; \hfill\break\indent spectral method} \begin{abstract} In this article we consider an energy balance climate model. For a given ice line, we use spectral method to derive an approximation of the solution. Then we propose a method to update the ice line and to derive an updated approximation of the solution. We compare the difference between the approximation with fixed ice line and the approximation with updated ice line by looking at the temperature profile at some specific locations and times. The significance of the method to update the ice line is that it is model free. Therefore, it can be used in other climate models. \end{abstract} \maketitle \numberwithin{equation}{section} \allowdisplaybreaks \section{Introduction} Zonal energy balance climate models(EBMs) were introduced by Russian climatologist Budyko \cite{MIBO} and Sellers of the University of Arizona \cite{WDSS} in the 1960s independently and have been studied extensively by North and his coworkers in the 1970s and 1980s \cite{EBM21,GRN,GRN1,GRN2,EBM7,GNJC} and other scientists and mathematicians recently \cite{FDHG,EBM9,EBM52,EBM20,KKT1,JWRM,JWEW,EBM16}. In this paper we adopt a version of zonal EBM from \cite{HKHE} (see also \cite{EBM21,FDHG,EBM9,EBM52,GRN,GRN1,GRN2,EBM20,GNJC}). For simplicity, we assume that all functions depending on latitude are symmetric about equator. Thus we need only to consider the Northern Hemisphere $(0<\theta<\frac{\pi}{2})$. Assume the spherical climate system varies with latitude but is uniform along lines of constant latitude. At any time $t$, let $T(t,x)$ be the mean annual surface temperature($^oC$) on the circle of latitude $\theta$ with $x=\sin\theta$. Let $C$ be the effective heat capacity of the system; that is, the energy needed to raise the temperature of Earth's surface by one degree Celsius per square meter. We consider the following diffusive zonally symmetric mean annual energy balance climate model (see also \cite{EBM21,FDHG,EBM9,HKHE,GRN,GRN1,GRN2,EBM20,GNJC}): \begin{equation}\label{eq:01} C\frac{\partial T}{\partial t} =D\frac{\partial }{\partial x}(1-x^2)\frac{\partial T}{\partial x} +QS(x)\beta(x,\mu)-(A+BT) \end{equation} where $Q$ is the globally averaged solar incident flux, $D$ is the horizontal heat diffusion coefficient and $S(x)$ is the normalized mean annual distribution of solar radiation. $I(T)=A+BT$ is the outgoing infrared radiation flux with $A$ and $B$ being empirical radiation coefficients. $\beta(x,\mu)=1-\alpha(x, \mu)$ is the co-albedo with $\alpha(x,\mu)$ being the planetary albedo, where $\mu$ is the sine of the latitude of the ice line. To solve \eqref{eq:01}, because of the symmetric assumption, it is reasonable to couple it with the following homogeneous boundary conditions at $x=0$ and $x=1$ as done in \cite{FDHG,EBM9,GRN,GRN1}, etc. \begin{equation}\label{eq:02} (1-x^2)^{1/2}\frac{\partial T}{\partial x}\big|_{x=0,1}=0 \end{equation} which means there is no heat transport across the equator and at the pole. We can see that several parameters are involved in the model. The variation of these prameters will change the temperature profile. But the change of these parameters is due to the physical changes on Earth, in the atmosphere, and other solar variations. It is believed and partially supported by geologic evidence \cite{EBM16} that \begin{quote} the solar variations of the past few million years are too small to account for the large fluctuations in climate such as glaciation/interglaciation events. Therefore, climate feedbacks, like those due to greenhouse gases and ice albedo feedback, provide explanations for these major changes and understanding the transient climate response to the feedbacks is of central concern in climate science. \end{quote} Here we are concern with the effects of ice line on the temperature profile. Since $\beta(x, \mu)$ depends on the location of ice line, as the solution of \eqref{eq:01}, Earth's temperature will depend on the location of the ice line as well. But the ice line is not fixed. Indeed, the higher the global temperature is, the less the ice cover is likely to extend, and thus the lower the albedo. A lower albedo leads in turn to a higher global temperature, and so on and so forth until all the ice is melted. In this paper we are going to introduce a method to update the location of the ice line and to investigate the effect of the movement of the ice line on the temperature profile. \section{A method for updating the location of iceline} In some primary energy balance climate models, it is assumed that Earth's surface is either ice covered or ice free and there is only one ice line (e.g. see \cite{MIBO,EBM21,FDHG,RMCL,GRN,GRN1,EBM7, JWEW,EBM16}, etc.). It is well-known that ice and snow have much higher albedos than bare soil or water. Therefore, the location of the ice line will determine how much the solar radiation will be absorbed by Earth. This in turn will affect Earth's temperature profile. But we know that the location of the ice line is changing with time. How to reflect the effect of the ice line change on Earth's temperature profile is an important question to answer. In \cite{JWEW,EBM16}, the authors proposed a method to investigate the movement of the ice line in the framework of manifold. In this paper we propose a method to update the ice line. To investigate the effect of the movement of the ice line on the temperature profile, we first set a reasonable value of $\mu$ in \eqref{eq:01}, then we use spectral method, that is, to represent the exact solution by a rapidly convergent sequence of approximations with each succeeding term incorporating information from successively smaller scales \cite{GRN1}, to solve \eqref{eq:01} as follows: Let $P_n(x)$ be Legendre polynomial of order $n$, then $\{P_n(x), n\;\text{even}\;\}$ forms an orthogonal basis for functions defined on $(0,\;1)$ and satisfies boundary condition \eqref{eq:02}. Assume that \begin{equation}\label{eq:21} T(t,x)=\sum_{n=0}^{\infty}T_n(t)P_n(x),\quad n \text{ even}, \end{equation} Then, $T(t,x)$ satisfies boundary condition \eqref{eq:02} term by term and \[ \frac{\partial T}{\partial t}(t,x)=\sum_{n=0}^{\infty}\frac{dT_n}{dt}(t)P_n(x). \] Since each term in the expansion \eqref{eq:21} of $T(t,x)$ satisfies boundary condition \eqref{eq:02}, when we truncate \eqref{eq:21} for any given even integer, \eqref{eq:02} is always satisfied. Insert these expansions into \eqref{eq:01}; multiply by $P_n(x)$ and integrate with respect to $x$ from $0$ to $1$. Making use of the orthogonality of $P_n(x)$ we obtain \begin{equation}\label{eq:12} C\frac{dT_n}{dt}+\left(n(n+1)D+B\right)T_n+A\delta_{n0}=QH_n, \end{equation} where \[ H_n(\mu)=(2n+1)\int_0^1S(x)\beta(x,\mu)P_n(x)dx. \] To solve \eqref{eq:12}, we need initial condition for each $T_n(t)$. Let $f(x)$ be the initial temperature distribution. Assume that $f(x)$ can be written as \[ f(x)=\sum_{n=0}^{\infty}f_nP_n(x), \] where \[ f_n=(2n+1)\int_0^1f(x)P_n(x)dx, \] then from \[ T(0,x)=\sum_{n=0}^{\infty}T_n(0)P_n(x)=f(x)=\sum_{n=0}^{\infty}f_nP_n(x), \] we have $T_n(0)=f_n$. For any given positive even integer $N$, we solve \eqref{eq:12} for $n=0, 2,4, \cdots, N$ to find $T_n$, then we find an N-mode approximation of $T(t,x)$ given by \[ T^N(t,x)\approx \sum_{n=0}^{N}T_n(t)P_n(x),\quad n \text{ even}. \] Assume that Earth's surface is either water (ice free) or ice covered and that there is only one ice line, $\mu$. Let \[ T(\mu)=T_c=\rho^o C \] be the critical temperature for the formation of ice. That is, \begin{gather*} T>\rho,\text{ no ice present,}\\ T<\rho,\text{ ice present.} \end{gather*} Let $M$ be a given positive integer and we solve \begin{equation}\label{eq:03} T^N(M,x)=\rho \end{equation} for $x$, say, $x=x_0$. Then we take $\mu=x_0$ in \eqref{eq:01} and repeat above process. Continue this process as many times as you like we will see the effect of the movement of the ice line on the temperature profile. \section{Execution of the method} As an example to demonstrate the method proposed in Section 2, following \cite{GRN2}, we take $A=203.3 W/m^2$, $B=2.09 W/m^{2\;o}C$, and $D=0.6487$. Following \cite{GKJL}, we take $Q=340.5$. It is known that the global heat capacity of Earth is no less than $2.08\times 10^8$ $J/m^{2o}C$. For demonstration purpose, we take $C=2.08\times 10^8$ $J/m^{2o}C$. Following \cite{GRN,GRN1,EBM7,EBM16}, we take $S(x)=1-0.482P_2(x)$ with \[ P_2(x)=\frac{1}{2}(3x^2-1). \] For simplicity we take $N=2$. Since the coefficients $T_n$ fall off very rapidly due to the factor $n(n+1)D+B$ in \eqref{eq:12}, we expect that a good approximation can be obtained from only the first two terms. For $n=0$, \eqref{eq:12} becomes \begin{equation}\label{eq:13} C\frac{dT_0}{dt}+BT_0=\rho, \end{equation} where $\rho=QH_0-A$ with $H_0$ being \begin{equation}\label{eq:013} H_0(\mu)=\int_0^1S(x)\beta(x,\mu)dx. \end{equation} Its general solution is \begin{equation}\label{eq:03b} T_0(t)=\frac{\rho}{B}+K_0e^{-\frac{B}{C}t}, \end{equation} where $K_0$ is an arbitrary constant to be determined from initial condition. For $n=2$, \eqref{eq:12} becomes \begin{equation}\label{eq:14} C\frac{dT_2}{dt}+(6D+B)T_2=QH_2, \end{equation} where \[H_2(\mu)=5\int_0^1S(x)\beta(x,\mu)P_2(x)dx.\] Its general solution is \begin{equation}\label{eq:04} T_2(t)=\frac{QH_2}{6D+B}+K_2e^{-\frac{6D+B}{C}t}, \end{equation} where $K_2$ is another arbitrary constant to be determined. After we obtain $T_0(t)$ and $T_2(t)$, we have a two-mode approximation of $T(t,x)$: \begin{equation}\label{eq:05} T^2(t,x)\approx T_0(t)+T_2(t)P_2(x), \end{equation} where $T_0(t)$ and $T_2(t)$ are given by \eqref{eq:03b} and \eqref{eq:04}. We take initial location of the ice line to be $\mu=0.95$ (same as in \cite{GRN1}). Same as in \cite{EBM21,GRN,GRN1,GRN2,EBM7,EBM20}, we adopt an ice cap parameterization that is due to Budyko \cite{MIBO}, namely, we take the critical temperature to be $-10^oC$. As for the co-albedo $\beta(x,\mu)$ we take it as a step function as done in \cite{MIBO,GRN,EBM7}, etc. \[ \beta(x,\mu)= \begin{cases} 0.38, & x>\mu,\\ 0.68, & x<\mu. \end{cases} \] Then by direct computations we have \begin{gather*} H_0=0.671697,\quad \rho=QH_0-A=25.412829\,,\\ T_0(t)=12.159248+K_0e^{-1.004808\times 10^{-8}t}. \end{gather*} To determine the constant $K_0$ we need an initial condition. From \cite{WBST} we know that the global mean annual temperature average for 2000s is $14.51^oC$. It is reasonable to take it as $T_0(0)$ since the $T_0$ in the expansion of the equilibrium solution is the global average temperature. By taking it as $T_0(0)$ we find that $K_0=2.350752$. Then \[ T_0(t)=12.159248+2.350752e^{-1.004808\times 10^{-8}t}. \] For the given values of the parameters, by direct computations, we have \begin{gather*} H_2=-0.366150,\\ T_2(t)=-20.8408+K_2e^{-2.876058\times 10^{-8}t}. \end{gather*} Again we need to determine the constant $K_2$. We take $T_2(0)$ to be the $T_2$ corresponding the equilibrium solution. From \cite{GRN2} we have $T_2=-28^oC$. Then, $K_2=-7.1592$. Therefore, \[ T_2(t)=-20.8408-7.1592e^{-2.876058\times 10^{-8}t}. \] and the two-mode approximation of $T(t,x)$ is \begin{equation}\label{eq:t201} \begin{split} \widehat{T}^2(t,x) &=T_0(t)+T_2(t)P_2(x)\\ &=12.159248+2.350752e^{-1.004808\times 10^{-8}t}\\ &\quad +(-20.8408-7.1592e^{-2.876058\times 10^{-8}t})\frac{1}{2}(3x^2-1)\\ &=22.579648+2.350752e^{-1.004808\times 10^{-8}t}+3.5796e^{-2.876058\times 10^{-8}t}\\ &\quad -31.2612x^2-10.7388x^2e^{-2.876058\times 10^{-8}t}. \end{split} \end{equation} Now we take $M=1000$ and solve \[ \widehat{T}^2(1000,\mu)=-10, \] we obtain $\mu=0.957553$. Now we can see that the ice line is moving up a little. Then we take the new $\mu$ in \eqref{eq:013} and obtain new \[ T_0(t)=12.3745+K_0e^{-1.004808\times 10^{-8}t}. \] This time, to determine the constant $K_0$, we take $\widehat{T}^2(1000,x)$ as initial condition. Then by direct computation, we found $T_0(0)=14.51$ and $K_0=2.1355$. Therefore, we have \[ T_0(t)=12.3745+2.1355e^{-1.004808\times 10^{-8}t}. \] Continuing our computations, we obtain \[ T_2(t)=-20.5157+K_2e^{-2.876058\times 10^{-8}t}. \] From the initial condition, we have $T_2(0)=-27.9998$. Then $K_2=-7.4841$ and \[ T_2(t)=-20.5157-7.4841e^{-2.876058\times 10^{-8}t}. \] Therefore, the new two-mode approximation of $T(t,x)$ is \begin{equation}\label{eq:t203} \begin{split} \widetilde{T}^2(t,x) &=T_0(t)+T_2(t)P_2(x)\\ &=12.3745+2.1355e^{-1.004808\times 10^{-8}t}\\ &\quad +(-20.5157-7.4841e^{-2.876058\times 10^{-8}t})\frac{1}{2}(3x^2-1)\\ &=22.63235+2.1355e^{-1.004808\times 10^{-8}t}+3.74205e^{-2.876058\times 10^{-8}t}\\ &\quad-30.77355x^2-11.22615x^2e^{-2.876058\times 10^{-8}t}. \end{split} \end{equation} Continuing to update the ice line. This time we take $M=10^5$ and solving \[\widetilde{T}^2(10^5,\mu)=-10\] gives us a new $\mu=0.957761$. By repeating above process, we obtain a new two-mode approximation of $T(t,x)$, \begin{equation}\label{eq:t202} \begin{split} \Bar{T}^2(t,x) &=T_0(t)+T_2(t)P_2(x)=12.3804+2.1275e^{-1.004808\times 10^{-8}t}\\ &\quad +(-20.5067-7.4716e^{-2.876058\times 10^{-8}t})\frac{1}{2}(3x^2-1)\\ &=22.63375+2.1275e^{-1.004808\times 10^{-8}t}+3.7358e^{-2.876058\times 10^{-8}t}\\ &\quad -30.76005x^2-11.2074x^2e^{-2.876058\times 10^{-8}t}. \end{split} \end{equation} This process can be continued. After we got this we can write out an updated approximation of $T(t,x)$ as a piecewise defined function as follows: \begin{equation}\label{eq:t204} T^2(t, x)= \begin{cases} 22.579648+2.350752e^{-1.004808\times 10^{-8}t}+3.5796e^{-2.876058\times 10^{-8}t}\\ -31.2612x^2-10.7388x^2e^{-2.876058\times 10^{-8}t}, \\ \quad \text{if }0\leq t\leq 10^3,\; 0(10^3+10^5),\; 0( 10^3+10^5). \end{cases} \] To compare the difference between the updated approximation \eqref{eq:t204} and the original approximation \eqref{eq:t201}, we still take New York City as an example. We take time $t=10^8$, by using \eqref{eq:t201}, we obtain \[ \widehat T^2(10^8, NYC)=10.081652^oC, \] and by using \eqref{eq:t204}, we obtain \[ T^2(10^8,NYC)=10.265468^oC. \] We can see that $ T^2(10^8,NYC)$ is a little higher than $\widehat T^2(10^8,NYC)$. We believe that $ T^2(10^8,NYC)$ should be a better estimate since it used the updated ice line. In general, for any location $x$, at time $t=10^8$, using \eqref{eq:t201} we obtain \[ \widehat T^2(10^8,x)=23.642027-31.866401x^2. \] Using \eqref{eq:t204}, the temperature profile is given by \[ T^2(10^{8}, x)=23.624605-31.393497x^2. \] We can see that there are some differences between these two temperature profiles although they are not very much. It is interesting to notice that, at $x=0$, that is, at the equator, \[ \widehat T^2(10^8,0)=23.642027,\quad T^2(10^{8}, 0)=23.624605. \] We can see that $\widehat T^2(10^8,0)>T^2(10^{8}, 0)$. But, at $x=1$, that is, at North pole, \[ \widehat T^2(10^8,1)=-8.224374,\;\;T^2(10^{8}, 1)=-7.768892. \] We can see that $\widehat T^2(10^8,1)