Cumulative Departure Model of the Cryosphere During the Pleistocene

A mathematical model is developed to describe changes in ice volume in the cryosphere. Modeling the cryosphere may be useful in assessing future climate impacts currently captured by global circulation models (GCMs) by providing an opportunity to validate GCMs. Leveraging the dominating effects of freezing and thawing in the cryosphere to simplify relevant heat transport equations allows for the derivation of a mathematical model that can be solved exactly. Such exact solutions are useful in investigating other climatic components that may be similarly analyzed for possible GCM validation. The current trend in GCM advancement is to increase the complexity and sophistication of the various heat transport effects that are represented in the governing mathematical model in cumulative form as the heat forcing function. In this paper, simplified models are developed whose solution can be directly compared with available data forms representing temperature and ice volume during the Pleistocene. With careful integration of the Pleistocene temperature term in the mathematical solution, the well-known cumulative departure method can be resolved from the mathematical solution using a two-term expansion of the corresponding Taylor series. This simplification is shown to be a good approximation of the Pleistocene ice volume for given Pleistocene temperatures. DOI: 10.1061/(ASCE)CR.1943-5495 .0000071. This work is made available under the terms of the Creative Commons Attribution 4.0 International license, http://creativecommons.org/licenses/by/4.0/. Author keywords: Global climate change; Cryosphere; Mathematical model; Phase change; Differential equation model; Cumulative departure method. Professor, Dept. of Mathematical Sciences, United States Military Academy, 626 Swift Rd., West Point, NY 10996 (corresponding author). E-mail: tedhromadka@yahoo.com Associate Professor, Dept. of Mathematical Sciences, United States Military Academy, 626 Swift Rd., West Point, NY 10996. E-mail: Doug.McInvale@usma.edu Instructor, Dept. of Mathematical Sciences, United States Military Academy, 626 Swift Rd., West Point, NY 10996. E-mail: Benjamin.Gatzke@usma.edu Department Head, Dept. of Mathematical Sciences, United States Military Academy, 626 Swift Rd., West Point, NY 10996. E-mail: Michael.Phillips@usma.edu Scientist, Hromadka and Associates, 29809 Santa Margarita Parkway, Suite 102, Rancho Santa Margarita, CA 92688. E-mail: bespinosa@sbcglobal.net Note. This manuscript was submitted on February 14, 2013; approved on March 12, 2014; published online on April 30, 2014. Discussion period open until September 30, 2014; separate discussions must be submitted for individual papers. This technical note is part of the Journal of Cold Regions Engineering, © ASCE, ISSN 0887-381X/06014002(14)/$25.00. © ASCE 06014002-1 J. Cold Reg. Eng. Journal of Cold Regions Engineering D ow nl oa de d fr om a sc el ib ra ry .o rg b y 98 .1 89 .2 14 .2 30 o n 04 /3 0/ 14 . C op yr ig ht A SC E . F or p er so na l u se o nl y; a ll ri gh ts r es er ve d.


Introduction
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).

Formulation of Mathematical Model
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. (1993Guymon et al. ( , 1981aGuymon et al. ( , b, 1980Guymon et al. ( , 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. 1999Petit et al. , 2001;;Jouzel et al. 1987Jouzel et al. , 1993Jouzel et al. , 1996Jouzel et al. , 2007a, b) , b) and a ratio of oxygen isotopes given by d18O (see Petit et al. 1999Petit et al. , 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. (1999Petit et al. ( , 2001)), Jouzel et al. (1987Jouzel et al. ( , 1993Jouzel et al. ( , 1996Jouzel et al. ( , 2007a, b), b)] to estimates of ice volume [as derived from proxy data, such as the ratio d18O; see Petit et al. (1999Petit et al. ( , 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 Proxy for ice volume: 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.

Data Types and Use of Proxy Relationships
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. 1999Petit et al. , 2001;;Jouzel et al. 1987Jouzel et al. , 1993Jouzel et al. , 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. 1999Petit et al. , 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. (1984Imbrie et al. ( , 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.

Mathematical Analog
In Hromadka [2012; Eq. ( 7)], a mathematical analog of changes in volumetric ice in the cryosphere is developed given by where 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 , gives where 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 cal=g, which suggests that 0 < r < 80 cal=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 0 = parameter cluster, Letting u ¼ 1 − kiðtÞ, then du ¼ −ki 0 ðtÞdt and Integrating from model time 0 to t gives Define GðtÞ by where 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) gives where ið0Þ ¼ 1 as the initial condition and Gð0Þ is defined to be zero.Multiplying by (k=r 0 ) and exponentiating Cumulative Departure Method Assume that heat contributed towards phase change, hðtÞ, is a linear function of a relevant temperature data source, hðtÞ ¼ aT þ b, with fa; bg = constants, then the mean value of hðtÞ over the simulation is a function of T avg , and where h avg ¼ aT avg þ b, and hðtÞ − h avg is the "departure" (from T 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Þ ¼ ∫ t s¼0 ½hðsÞds − h avg t or where, from Eq. ( 7), GðtÞ ¼ ∫ t 0 hðsÞds and so combining Eqs. ( 10) and ( 13): 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., aT avg þ b ≈ 0) and Eq. ( 15) simplifies to where 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 0 > 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, gives where D = cluster of parameters (i.e., example parameters) given by where r 0 ¼ rI 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.

Modeling Demonstrations
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 longterm 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.

Discussion
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   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.

Conclusions
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.