Cumulative Departure Model of the Cryosphere During the Pleistocene

The recent literature regarding assessment of the quality of information produced from global climate models [or global circulation models (GCMs)] has been mixed as far as the successes achieved by use of GCM modeling (Randall et al. 2007; Anagnostopoulos et al. 2010; Bounoua et al. 2010). Notwithstanding these mixed assessments, use of GCMs appears to be a main-stream course in the analysis of climate change effects to develop public policy. Because the direction of GCM development is to include additional components of the heat budget and its several transport mechanisms, such advancement may need to address the well-known uncertainty problems that have been identified in the focused modeling efforts of those individual transport mechanisms (e.g., Hromadka 2012). However, some aspects of the climate may be more predictable, depending on the time scale involved and the primary heat budget elements involved. For example, a model of the heat budget associated with the cryosphere may be more tractable because of the dominating effects of the latent heat of phase change in freezing and thawing (e.g., Dyurgerov and Meier 1999; Bounoua et al. 2010; Bamber and Payne 2004). In the current work, a mathematical model of ice volume in the cryosphere is developed and a closed-form solution is obtained and subsequently numerically modeled using the well-known cumulative departure method (CDM).

Similar to the underpinnings of the Stefan model (e.g., Alexiades and Solomon 1993) for freezing and thawing soils or other medium, the heat budget relationship can be simplified because of the dominating effects of freezing and thawing of water in the cryosphere. For example, in the well-known Stefan solution (Alexiades and Solomon 1993) for estimating the freezing front depth in a column, steady-state conditions are assumed above and below the freezing front, with heat exchange being dominated by the phase change process at the freezing front. In Guymon et al. (1993, 1981a, b, 1980, 1984), Guymon and Hromadka (1980), Hromadka (1986a, b, c, d), and Hromadka and Guymon (1982, 1984, 1985), finite-element and complex variable boundary element method computer models of freezing and thawing soils in algid climates were developed for two- and three-dimensional problems, including impacts from groundwater flow regimes and with isothermal phase change in moist soils. These models are based on the dominating effects of the phase change process in the heat transport budget, and have been applied to problems involving soil water phase change and water-ice interface phase change effects in long-term simulations.

Data forms available that may be used to describe heat budget characteristics include isotopes of oxygen and hydrogen. Commonly used ratios of isotopes are used to provide proxy information of air temperature and ice volume, such as a ratio of hydrogen isotopes given by $dD$ (see Petit et al. 1999, 2001; Jouzel et al. 1987, 1993, 1996, 2007a, b) and a ratio of oxygen isotopes given by $d18O$ (see Petit et al. 1999, 2001; Lisiecki and Raymo 2005a). Equations for these proxy data types are listed in Table 1, with terms and relationships defined in Table 2:

The proposed conceptual model will directly relate estimates of temperature found in the literature [as derived from proxy data, such as the ratio $dD$; for example, see Petit et al. (1999, 2001), Jouzel et al. (1987, 1993, 1996, 2007a, b)] to estimates of ice volume [as derived from proxy data, such as the ratio $d18O$; see Petit et al. (1999, 2001), Lisiecki and Raymo (2005b)]. Given a mathematical analog of temperature and ice volume in differential equation form, an exact solution is obtained. For situations not suitable for direct solution, a central-difference finite-difference method analog, using computer program EXCEL, may be developed (e.g., Hromadka 2012).

Model calibration is made to normalized data of frozen material, $I(t)$, for the Pleistocene time period. To simplify the calibration process, the governing differential equation is normalized to reduce the number of parameters. Some modeling parameters are directly calibrated to the selected proxy data representing the cumulative heat magnitude function, $H(t)$. The other modeling parameters are calibrated to fit the modeling outcomes to the selected measure of frozen material $I(t)$.

The heat transport effects that impact the problem domain include solar and internal heat within the planet, circulation of the atmosphere and oceans with corresponding interface, and convective heat transport effects, among others. All of these components interface and integrate together and are represented by the term [Heat]. [The notation [Heat] includes all the components of heat transport that effect phase change in the cryosphere. The total heat budget can be found in Bitz and Marshall (2012), among others]. [Heat] includes all the thermal balance elements of the heat budget, with magnitude, $H(t)$. The function $H(t)$ is determined by calibration to the historic temperature $T(t)$ function developed from historic data. In the current approach, $H(t)$ is developed by examining what $H(t)$ did in the past with respect the $T(t)$ function. This modeling approach of building $H(t)$ based on what $H(t)$ did in the past may sidestep many complexities that are involved in assembling approximations of the various subprocesses and feedbacks that form the climate and associated boundary conditions. Furthermore, such an $H(t)$ can still be modified to represent changes in the environment. Inferences can then be drawn that relate to the climate outside of the cryosphere, given the effects within the cryosphere. The resulting model framework can be used as model verification for GCMs, under suitable boundary and initial conditions, for the effects on the cryosphere, such as the Arctic and Antarctic regions.

The solution to the initial basic conceptual model relates an input forcing function, $H(t)$, to the output $I(t)$ function. The input function is assumed to follow the trends of reported data, such as ice core temperature estimates (Jouzel et al. 2007a; Kawamura et al. 2007b; Petit et al. 2001). Under assumed “natural conditions,” the output function, $I(t)$, is assumed to also follow the data trends that relate frozen material measure to proxy data such as $d18O$.

Because normalized data are used, the forcing function, $H(t)$, is linearly related to the selected temperature data proxy for the time domain of interest. That is, for some constants $a_1$ and $a_2$, and assuming $H(t)=a_1+a_{2}T(t)$, then when transformed into normalized $N(0,1)$ form, both functions $H(t)$ and $T(t)$ plot identically. Consequently, one can work directly with the $N(0,1)$ normalized form of a selected $T(t)$ realization. Furthermore, the estimated temperature $T(t)$ is typically given as a linear function of the proxy data, $dD$ (e.g., Petit et al. 1999, 2001; Jouzel et al. 1987, 1993, 1996). Consequently, the model formulation necessarily is using a linear transform of $dD$ as the forcing function, and therefore the $N(0,1)$ transform of the available $dD$ source data serves to describe the forcing function.

Some papers suggest that ice volume estimates may be developed as a linear function of $d18O$ (e.g., Petit et al. 1999, 2001; Zachos et al. 2001), and other papers suggest that ice volume estimates may be developed using a nonlinear function of $d18O$ (e.g., Mix and Ruddiman 1984). Assuming that the relationship is linear for the range of conditions contemplated in the modeling formulation, the $N(0,1)$ transform of the available $d18O$ source data describes the ice volume trends during the Pleistocene time period. (If the relationship is nonlinear, then the $N(0,1)$ transform of the $d18O$ source data may require further considerations, such as including other statistical moments or using another choice of transformation.) In the current paper, $d18O$ data (Kawamura et al. 2007a; Petit et al. 2001) are used as a proxy for representing the trends of ice volume, which becomes the target for the calibration of the conceptual model formulation.

From the preceding discussion, it is seen that the general solution to the conceptual model operates on estimates of the magnitude of heat, $H(t)$, assumed to be a linear function of temperature. In turn, these are assumed to be a linear function of $dD$ and produce an output of ice measure that is calibrated to estimates of ice volume. Ice volume is assumed to be a linear function of $d18O$ (as stated previously, the preceding development can consider a nonlinear transformation of $d18O$, or other proxy data). Hence, the model solution is a mapping (i.e., a relationship between the two measured variables considered) from $dD$ to $d180$.

The proxy data types considered in this paper include those listed in Table 1. Each type of proxy shown in bold in Table 2 were collected and processed by partitioning the study time domain into 1,000-year intervals [the selection of a 1000-year partition is consistent with the partition size used in Imbrie et al. (1984, 2011) and Imbrie and Imbrie (1980) in a study of phase-space modeling for Pleistocene ice volume]. The resulting partition was modified to better fit the source data relative to the maximum and minimum data to capture the respective peak and minimal values. The modified partitioned data were then transformed using the standard normal $N(0,1)$ transformation. In such a normalized form, calibration of the forcing function $H(t)$ to a selected proxy and calibration of the governing PDE solution $I(t)$ is made simpler, as fewer parameters are involved. It is observed in Hromadka (2012) that the different data sources for a particular variable, such as estimated temperature or estimated ice volume, do not agree well. Consequently, any model using such a data source as input (or output for calibration) would be impacted by such differences. Assuming that the assemblage of identical data types are equally likely realizations of a stochastic process, then the forcing function may be modeled as a set of input realizations (i.e., estimated temperature), each equally likely, and the model then calibrated to each of the assembled output realizations (i.e., estimated ice volume or measure), each equally likely.

In Hromadka [2012; Eq. (7)], a mathematical analog of changes in volumetric ice in the cryosphere is developed given bywhere $H(t)$ = total heat into the system affecting phase change in the cryosphere; $I(t)$ = total volumetric ice in the cryosphere including glaciers, tundra ice, snow, and other forms of water subject to freeze/thaw effects; $I_0$ = an initial condition of $I(t)$, contemplated as a local maximum value such as during a glacial period; $k$ = lumped parameter representing that portion of heat returned due to phase change effects, such as reflected heat due to snow, among other factors, where $0<k<1$; $r$ = lumped parameter representing the area-averaged latent heat of fusion as averaged throughout the cryosphere; and $t$ = model time.

For details regarding the underpinnings of Eq. (1), the reader is referred to Hromadka (2012).

The ratio $I(t)/I_0$ in Eq. (1) represents the impact of changing geometry during phase change effects, such as growing glaciers, or changing boundary lengths with respect to interface with ocean currents, among other effects. In Hromadka (2012), $i(t)=I(t)/I_0$ is used with an exponent, $p$, such that $i{(t)}^p$ is considered in Eq. (1). For the current work, however, $p=1$ is retained.

Dividing $dI(t)/dt$ by model initial condition, $I(t=0)=I_0$, giveswhere the initial condition for $I(t=0)=I_0$, or $i(t=0)=1$. The literature (Ruddiman 2001) indicates that the minimum value of $i(t)$ during the Pleistocene period may have a value of about 0.30. For pure water, the latent heat of fusion has value $80\mathrm{cal}/\mathrm{g}$, which suggests that $0<r<80\mathrm{cal}/\mathrm{g}$ may be an appropriate range for the parameter, $r$.

It is useful to divide Eq. (2) by a reference value for heat given by $H_0$. Then letting $h(t)=H(t)/H_0$where $r^{\prime }$ = parameter cluster, $r^{\prime }=(r{I}_0/H_0)$.

Rearranging terms

Letting $u=1-ki(t)$, then $du=-k{i}^{\prime }(t)dt$ and

Integrating from model time 0 to $t$ gives

Define $G(t)$ bywhere $G(t)$ = total heat applied to the cryosphere between model time zero and $t$; and where $G_0=G(s=0)$ = initial condition value.

Then combining Eqs. (6) and (7) giveswhere $i(0)=1$ as the initial condition and $G(0)$ is defined to be zero. Multiplying by ($k/r^{\prime }$) and exponentiatingor

Assume that heat contributed towards phase change, $h(t)$, is a linear function of a relevant temperature data source, $h(t)=aT+b$, with $\{a,b\}$ = constants, then the mean value of $h(t)$ over the simulation is a function of $T_{\mathrm{avg}}$, andwhere $h_{\mathrm{avg}}=a{T}_{\mathrm{avg}}+b$, and $h(t)-h_{\mathrm{avg}}$ is the “departure” (from $T_{\mathrm{avg}}$) of heat, evaluated at time, $t$. The cumulative departure of heat is given by $C_h(t)$, where from Eq. (11)where $C_h(t)$ = associated cumulative departure of the reference temperature.

From Eq. (12), $C_h(t)={\int }_{s=0}^t[h(s)ds]-h_{\mathrm{avg}}t$ orwhere, from Eq. (7), $G(t)={\int }_0^{t}h(s)ds$ and so combining Eqs. (10) and (13):or

In Eq. (15), the second exponential term indicates a long-term trend of aggradation or degradation of $i(t)$. If there is no evidence of such long-term aggradation or degradation of $i(t)$, then the argument of the second exponential term may be assumed zero (i.e., $a{T}_{\mathrm{avg}}+b\approx 0$) and Eq. (15) simplifies towhere $C_T(t)$ = outcome of the cumulative departure method applied to the reference temperature (it is noted that the reference temperature may be offset by an arbitrary constant and yet $C_T(t)$ is not changed in value).

Based on Eqs. (14) and (15) and the simplified version [Eq. (16)], the set of calibration parameters are {$k$, $r$, $a$, $b$} where $0<k<1$ and $r^{\prime }>0$.

Carrying forward the formulation of Eq. (16) [indicating no long-term aggradation or degradation of $i(t)$], the exponential function may be expanded into a Taylor series and neglecting terms of order 2 and higher, giveswhere $\mathfrak{D}$ = cluster of parameters (i.e., example parameters) given bywhere $r^{\prime }=r{I}_0/H_0$ [see text following Eq. (3)], and it is recalled that in this formulation long-term aggradation or degradation of $i(t)$ is assumed negligible.

In the analog of Eq. (19), the well-known cumulative departure method, as used in many fields of geoscience, is seen. As such, experience with using a formulation such as Eq. (19) will follow, as used in other areas of application of the CDM, such as groundwater budgeting, accounting, and other topics.

In this paper, two solution methods are developed for the mathematical model of cryosphere ice volume as given in Eq. (1), which is developed in detail in Hromadka (2012). First, an exact solution to the differential equation is developed and is shown in Eq. (10). Second, the exact solution is then resolved into a long-term background term and the variation about the long-term background term. This variation term is a variant of the CDM that is used to approximate other phenomenon in geoscience, such as rainfall trends or groundwater trends, among other topics. This background term is seen as the second exponential term of the full formulation given in Eq. (15), with the variation about the background term seen as the first exponential term in Eq. (10). By assuming that the Pleistocene ice volume trends seen most recently in the ice volume data records, such as the last 300,000 years or so, are not significantly influenced by a long-term ice volume aggradation or degradation, then the argument of the said second exponential term in Eq. (10) may be assumed to be approximately zero, reducing the formulation of Eq. (10) to the formulation of Eq. (16), which after some algebraic rearrangement of terms arrives at the final CDM formulation shown in Eqs. (17)–(20).

The simplified exact solution (normalized EXP) of Eq. (16) is compared with the finite-difference method (FDM) modeling results shown in the paper of Hromadka (2012) in Fig. 1. From Fig. 1, the two modeling approaches, namely, the one-term exponential version of the exact solution of Eq. (16) and the FDM compare closely. Shown in Fig. 2 are the results from the said two modeling approaches (shown in Fig. 1) plotted with the SPECtral MApping Project (SPECMAP) ice volume data for the same time period in the Pleistocene period of approximately the last 300,000 years. SPECMAP is NASA’s SPECtral MApping Project containing historical climatic time series data. In these modeling approaches, the same calibration parameter values are used as determined in Hromadka (2012).

In Fig. 3, the one-term version of the exact solution and the FDM results of Figs. 1 and 2 are compared with the results from the cumulative departure method (for the same parameter values). The three modeling approaches all compare closely. Fig. 4 compares the three modeling outcomes with the ice volume data.

All of these modeling approaches were developed using the computer program EXCEL. Also, all of these modeling exercises are suitable for undergraduate mathematical modeling application problems, such as those found at the U.S. Military Academy at West Point, New York.

The simplification of the exact solution [i.e., Eq. (16)] to the mathematical model of Eq. (1), resulting in the formulation seen in Eqs. (17)–(20), develops a new application of the CDM. The underpinnings of this particular application of the CDM is found in Eqs. (16)–(20) and involves a single parameter that is a cluster of the parameters used in the mathematical model of Eq. (1).

The modeling formulations all assume that the problem boundary condition at model time zero corresponds to a relative maximum value of ice volume where $i(t=0)=1$. Consequently, model ice volume values may exceed the value of 1 whenever the estimated ice volume exceeds the ice volume occurring at model time $t=0$. However, data show that volumetric ice volume does not fall to less than approximately 0.30 (e.g., Ruddiman 2001). Consequently, these bounds for the ice volume variable can be used to calibrate the single CDM parameter seen in Eq. (20).

The CDM model can be readily used to estimate the impacts of climate changes by providing the change in ice volume given an alteration in one of the modeling parameters seen in Eq. (20). For example, if a particular effect is concluded to have approximately a 5% increase in the heat returned [as seen in Eq. (1)], then the CDM parameter of Eq. (20) can be recalculated and used accordingly. The overall influence of the particular effect being considered may be estimated by “rewinding” the CDM model to some point in model time (perhaps some 200 years), just prior to the current time, and then evaluate the CDM model but with the parameter change implemented for only the subject investigated time period. Further research is needed to explore the possibilities available in using the CDM formulation in assessing some of the aspects of climate change.

A simplified model of ice volume in the cryosphere during the Pleistocene time period is developed using the well-known cumulative departure method. The starting point of the CDM model development is the differential equation formulation of cryosphere ice volume during the Pleistocene time period, as presented in Hromadka (2012). An exact solution is determined for this differential equation formulation, involving the product of two exponential terms, one of which represents the effects of significant long-term ice volume degradation or aggradation. Then, by assuming negligible long-term ice volume aggradation or degradation, the exact solution is simplified into a single-term exponential model that can be resolved into another application of the cumulative departure method, which is a modeling technique that is used in many other applications in geoscience. The resulting CDM formulation is a single parameter model based on a cluster of the several model parameters used in the original mathematical model formulation. Modeling results produced by the exact solution, and also the CDM, compare closely with the modeling results produced from a finite-difference method formulation given in Hromadka (2012). An approach to assessing possible climate change effects on the cryosphere ice volume is suggested. This approach adjusts the single CDM modeling parameter according to the parameter cluster underpinnings, and then rewinds the CDM model in model time and reevaluates using the revised parameter value. The various modeling approaches are tractable, assessable through standard computational programs, and has demonstrated utility for both practitioners and researchers.

Acknowledgements are paid to the United States Military Academy, West Point, New York, Department of Mathematical Sciences, for their support to the authors during this research. Also acknowledged are the several individuals who have participated in particular tasks in developing this paper, including but by no means limited to Rene Perez, Laura Hromadka, Michael Barton, and T.V. Hromadka III.