Towards a More Realistic Young-Earth Ice Sheet Model A Shallow, Isothermal Ice Dome with a Frozen Base

TOWARDS A MORE REALISTIC YOUNG-EARTH ICE SHEET MODEL A SHALLOW, ISOTHERMAL ICE DOME WITH A FROZEN BASE

Key words: ice core, ice dome, ice sheet, Larry Vardiman, M. W. Mahaffy, Ed Bueler, Ice Age, computer modeling

Introduction

In 1993 Larry Vardiman of the Institute for Creation Research published an analytical model for the rapid formation of thick ice sheets within the 4,500 years since the Genesis Flood (Vardiman, 1993, 1994, and 2001). Vardiman’s model assumed that post-Flood yearly accumulation decayed exponentially from a maximum value just after the Flood to today’s “slow and gradual” values. He implicitly assumed that at the center of the ice sheet the downward speed of the surface ice was proportional to the height of the ice sheet (Hebert, 2021, p. 179­). This assumption helped make the math tractable, but as Vardiman himself noted, this was a preliminary effort, and this assumption is not physically realistic. Hence there is a need for improved young-Earth ice sheet models, so as to make predictions and to compare and contrast those predictions with those made by uniformitarian models. In fact, glaciologists have already devised ice flow models which can be stripped of the usual uniformitarian assumptions and could conceivably be used for this purpose. In 1976 M. W. Mahaffy published a quasi-analytical ice sheet model which does not make the typical uniformitarian assumption that the ice sheet is in “steady state”, i.e. that it has a constant height. This paper combines Vardiman’s accumulation model with Mahaffy’s ice sheet model to obtain a more realistic estimate of the manner in which a 3,000-meter-thick ice sheet would grow over time. The results are then compared with results from Vardiman’s model. Limitations of the Mahaffy model are discussed, as are future avenues for creation ice core research. The next four sections provide an introduction to ice core basics. Readers already familiar with this information may wish to skip ahead to the section entitled “The Mahaffy Model”.

Ice Sheet Fundamentals

Glaciologists drill ice cores at geological features called ice divides. These features are named for the fact that ice on one side of the divide flows one way, and ice on the other side of the divide flows the other (Figure 1). For a perfectly flat surface, ice is deposited in horizontal layers. Assuming that no shifting of the ice divide during buildup has occurred, ice at the divide itself will flow neither to the right nor left, but will move straight down as additional layers are added. It is also possible for the ice to form a dome around the divide, so that ice flows radially outward from the divide, not just right and left. In fact, in this paper we will be modelling the formation of an ice dome.

            If we ignore possible isostatic adjustments, the underlying bedrock is locked in place. So a horizontal layer that reaches bedrock can no longer move downward. Furthermore, stresses on the ice layer will cause it to become thinner over time. If no melting occurs, the total number of annual layers will equal the number of years since the ice sheet began forming. Nevertheless, identification of annual layers is not trivial. Visibly distinguishable layers may not be present if annual accumulation is too small (Anonymous, 2020), as is the case in central Antarctica (Anonymous, 2021). And even when distinct bands are visible, multiple layers are deposited per year, and the number of these layers can vary from year to year (Alley and Koci, 1988; Alley, 1988). Hence one cannot naively assume each distinct band in an ice core represents an annual layer, as was inadvertently highlighted by creation critic Bill Nye in 2018 (Hebert, 2018a). Glaciologists make educated guesses as to how many bands should be grouped together and counted as a year (Alley et al., 1997, p. 26368), guesses that unfortunately are biased by circular reasoning and long-age assumptions (Oard, 2005; Hebert, 2014). Although uniformitarian glaciologists now think they can accurately count back tens of thousands of years with an accuracy of 10% or better (Meese et al., 1997, p. 26413, their Table 2), this was not always the case. Early pioneers (Hammer et al., 1978, p. 5, their Table I) in the field of ice core studies did not think it was possible to count back more than 200 years while maintaining an accuracy of 10%! This earlier, more pessimistic view was likely more realistic.  

            Ice sheet models usually ignore variations in density, which are very slight even for very thick ice sheets (Cuffey and Paterson, 2010, pp. 12-13). Hence, most ice sheet models assume incompressibility, that the ice has a constant density and volume (Paterson, 1980, p. 12; Cuffey and Paterson, 2010, p. 286).

            In the absence of melting or calving (ice breaking off a glacier’s edge), the volume of an ice layer remains constant over time. The layer becomes thinner, with a corresponding increase in the layer’s horizontal surface area, but in such a way that the total volume of the layer remains constant.

Stress Basics

Consider a small cube of ice with length, width, and height of dx, dy, and dz, and whose edges are aligned with the x, y, and z axes (Figure 2). Stresses are applied forces per unit area (units of Newtons per square meter, or Pascals) that tend to change a body’s shape and volume. For a cube of material, these stresses may be applied either perpendicular to or parallel to a cube face. Perpendicular stresses are called normal stresses, while parallel stresses are called shear stresses.

            The effect that a stress has upon an ice parcel depends on both the direction in which the force is pointing, as well as where the force is applied. For this reason, two indices are needed to specify a particular stress component , one to identify the direction in which the stress is pointing, and a second index to indicate on which pair of opposing cube faces the stress is being applied. The first subscript y in indicates a stress caused by a force pointing in the y-direction, and the x subscript indicates that the force is being applied to the two cube faces which are perpendicular to the x-axis, i.e., the front and rear faces of the cube in Figure 2.  is a shear stress. Likewise, the normal stress  is caused by a force pointing in the x-direction and applied to the faces that are perpendicular to the x-direction, i.e., the front and rear cube faces. As is demonstrated from these examples, represents a normal stress when , and a shear stress when .

            If the normal stress is the same (same size and direction) on two opposite cube faces, these stresses will tend to translate (move) the ice parcel through space, but not deform it. Likewise, if the same shear stresses are applied to both the top and bottom faces, the forces associated with those stresses will tend to translate the cube horizontally, but not deform it. In order for a stress to have a deformational effect on a small cube of ice, the stress must vary across the width of the cube.

            It is the difference in applied normal stresses that tends to compress or extend the parcel. It is fairly obvious that the difference in normal stresses between two opposite faces will tend to extend or compress the cube. Unbalanced shear stresses, on the other hand, tend to (Cuffey and Paterson, 2010, p. 680) deform the parcel (Figure 3).

            To calculate the manner in which an ice sheet deforms, it is necessary to sum up the stresses acting on a small ice parcel, much in the same way one sums up forces to find an object’s acceleration. For an incompressible fluid of density ρ, this requires the use of a partial differential equation called the Navier-Stokes equation (Bueler 2016, p. 2), which is the stress balance equivalent to Newton’s Second Law.

Because an ice sheet deforms so slowly, the acceleration (material time derivative of the momentum, or the time derivative of momentum for a moving parcel) in the Navier-Stokes equation will be zero, and forces will be balanced in all three directions. This is equivalent to solving the following three stress balance equations:

                                                                                  

Strain and Velocity

Normal strain is the fractional change in one of the linear dimensions of a material due to an applied normal stress, and the strain rate is the time derivative of the strain. If u, v, and w are the x, y, and z ice velocity components, it is fairly easy to show that the three normal strain rates (for small deformations) are given by (Cuffey and Paterson, 2010, p. 679)

                                                                                                                

Negative normal strains will compress a parcel, while positive normal strains will expand it. The shear strain rates (again, assuming the deformations are small) are given by (Paterson, 1980, p. 16; Cuffey and Paterson, 2010, p. 680):

                                                                                                                                     

Glen’s Flow Law

A very important relationship is Nye’s (1957) generalization of Glen’s experimentally-determined flow law (1955), which gives the relationship between stress and strain in ice, discussed below. Stresses can either change the volume of a body or its shape. Volumetric or hydrostatic stresses are associated with changes in volume, whereas deviatoric stresses or stress deviators are associated with changes in shape (but not volume). Since ice is essentially incompressible, the shape of an ice parcel can change, but not its volume. Hence one might expect that the flow law would be expressed in terms of deviatoric stresses, and that is indeed the case. The “stress deviator” is the part of the normal stress in excess of the mean stress. For an isotropic material, this mean stress is -P (although the hydrostatic pressure P is positive, the stress caused by P is -P, since hydrostatic stress is compressive). Hence, the “stress deviators” or “deviatoric stresses”are:

                                                                                                               

where δ_ij is the Kronecker delta, which is defined to be 1 when i = j, and 0 when i ≠ j. Because the Kronecker delta is equal to zero when i ≠ j, Equation implies that a shear stress is always deviatoric, which makes sense, since shear stresses change a parcel’s shape but not its volume. It is only the normal stresses (for which i = j and the Kronecker delta is 1) that can be expressed as the sum of volumetric and deviatoric stresses. As noted earlier, we want to express the flow law in terms of the deviatoric stresses. Since the ice is already assumed to be isotropic, the strain rates can be expected to be in the same direction and proportional to the stresses (Paterson, 1980, p. 15; Bueler, 2016, p. 6):

                                                                                                                             

where η is the ice viscosity. Analogous expressions hold for the other strain rates. Because the flow law is a physical property of the ice, its mathematical form should be independent of the choice of coordinate axes (Paterson, 1980, p. 15) . Nye (1957) proposed defining the effective strain rateand the effective shear stressas

                                                                                       

Since these quantities are the second invariants of the strain and stress tensors, respectively, they are unaffected by a rotation of the coordinate axes (Paterson, 1980, p. 16), and we can express Nye’s generalization of Glen’s Flow Law as

                                                                                                                                           

where A(T) is the temperature-dependent “softness” of the ice (warmer ice is softer than colder ice). The exponent n is normally taken to be 3 for glacial ice. Equations and imply that

                                                                                                                                                   

For n = 3, Equations and imply that  

                                                                                                                                             

Finally, this gives us the relationships

                                                                                     

with analogous relations holding for the other strain components.

Since the ice is assumed to have no rotational acceleration, the torques about any point in the ice must be zero. This condition is satisfied if the stress tensor is symmetric (Cuffey and Paterson, 2010, pp. 675-676)., so that By Equation , the strain rate tensor is also symmetric, with  

The Mahaffy Model

Mahaffy’s (1976) ice sheet model does not assume the ice to be in steady state. Hence, his model is of interest to creationists. Although Mahaffy considered the more general case of an ice sheet on non-level terrain, we will begin with the simpler case in which the bedrock under the ice is perfectly flat. We also treat the ice as isothermal, so that the ice softness A(T) is a constant A. Finally, we assume that there is no melting at the base of the ice sheet, so the ice is frozen to the bedrock.

Consider the ice sheet depicted in Figure 4. At a given time, the height h of the ice sheet can vary from place to place, as shown in the diagram. In other words, at a particular time, h is a function of the horizontal coordinates x and y. We will use H = h (0, 0) to denote the height at the center of the ice sheet, i.e. the divide location.

Consider also a thin column of ice with height h(x, y) and having a square base with sides of length 1 meter (Figure 5). We define q_x to be the volume of ice that flows horizontally (in the x-direction) into the column in a time dt.

                                                                                 

Note that if q(x) equals q(x+1), so that  = 0, there will be no net inflow or outflow of ice (in the x direction) into the column. On the other hand, if q(x+1) is less than q(x), so that , there will be a net inflow of ice, which will tend to cause the height of the ice column to increase. Likewise, a positive partial derivative corresponds to a net outflow of ice. Similar relationships hold for the y-direction. The net horizontal flow of ice into or out of the column will tend to change the height h, but so will the deposition of accumulating ice, indicated by . In general,  is a function of both position (x, y) and time t.

            Hence, the height h(x, y, t) of the ice is governed by the continuity equation:

                                                                                                                                 

Glaciologists have long assumed that the shear stress parallel to the bed is most important when discussing ice flow (Weertman, 1961, p. 954). Hence Mahaffey ignored deformation due to normal stresses, as well as deformation due to shear stresses in the vertical direction. In other words, he assumed that the ice sheet undergoes deformation due only to shear stresses in the x-y plane. This implies that the strain rates are:

                                                                                                 

Equation implies that the corresponding (deviatoric) shear and normal stresses are also zero. From Equation , it follows that . We also assume that vertical normal stress changes much more rapidly in the vertical direction than do the shear stresses in the horizontal directions:

                                                                                                                                 

The stress balance equations then reduce to

                                                                                                                               

Atmospheric pressure is negligible compared to the pressure inside the ice, so we may take the pressure at the surface to be P (z = h) = 0 without loss of generality. Integrating the last of Equation and applying this boundary condition give:

                                                                                                  

We further assume that the slope of the ice at the surface is small in both the x and y directions. This is a reasonable approximation at most locations in the ice, except at the edges, which can be quite steep. Also, as discussed later, normal stresses cannot really be neglected at the center of the ice sheet. Nye (1952, p. 86) argued that, for ice sheets with shallow surface slopes, the shear stresses depend mainly on the surface slope and the weight of the overlying ice. For shallow surface slopes, the shear stresses parallel to the ice surface are approximately horizontal, and we can approximate them as

                                                                                                                     

We have just implicitly made the “shallow ice approximation”. Equation and imply that

                                                                              

where

                                                                                                                                     

since all the other deviatoric stresses are zero.

Mahaffy also assumed that and . As can be seen from the flow lines in Figure 1, the vertical variations in the horizontal velocities u and v are much greater than the horizontal variation in the vertical velocity w, and this assumption is justified. After eliminating the partial derivatives of w with respect to x and y, we obtain

                                                                                                             

Inserting Equations into Equations , integrating twice with respect to z, and making use of the boundary conditions u(x, y, z = 0) = 0 and v(x, y, z = 0) = 0 yields the following expressions (Paterson, 1980, p. 43) for the horizontal velocity fluxes q_x and q_y:    

                                                                                               
where

                                                                                                                              

is the magnitude of the surface gradient. Since the greatest horizontal fluxes occur deep in the ice, A should be weighted in favor of the values of A deep in the ice. However, since we are treating the ice as isothermal, we do not need to worry about this complication, and we treat A as a constant. If we define

                                                                                                                                         

then we may combine Equations into a single expression for the flux vector :

                                                                                                                         

A Clever Solution

Bueler (2016) presents a very clever solution. Note that we can express the continuity equation in terms of the flux vector :

                                                                                                

Bueler (2016, p. 8) draws our attention to similarities between Equation and the 2-D heat equation:

                                                                                            

where f is a heat “source” term, ρ is density, c is specific heat, and κ is thermal conductivity. If we divide both sides of Equation by ρc, the result has the same mathematical form as Equation :

                                                                                                                               

where we have defined D and F as

                                                                                                                                                         

The choice of the letter D in Equation is appropriate, because the heat equation is actually a diffusion equation: any localized temperature ‘spike’ will spread out, with the peak becoming shorter and broader over time. Note that the same thing will be true for an ice peak; the ice flows in such a way that the ice peak becomes shorter and broader over time. We can thus treat the flow of the ice as a diffusion problem if we define our position and time-dependent ice “diffusivity” to be

                                                                                                                 

Bueler (2016) describes how Equation may be solved using an explicit, finite difference solution with adaptive time-stepping, and he provides a sample MATHLAB code for doing so. I “translated” his code into an IDL (Interactive Data Language) code I named mahaffy_vardiman_dome.pro, and used my version to simulate the growth of a large isothermal ice dome with a frozen base. Details and subtleties of the method are given in Bueler (2016).

Test Case

I applied Vardiman’s (1993, 1994, 2001) model for post-Flood ice deposition

                                                                                                                           

to an isothermal, symmetrical ice dome on perfectly flat bedrock. I chose parameter values relevant to the EPICA Dome C ice core. The accumulation rate is λ_H/τ = 2.84 centimeters of ice per year, the present-day ice accumulation at the ice core location (Parrenin et al., 2007, pp. 248). The parameter Φ, used to “scale up” the amount of Ice Age precipitation compared to the present-day value, was set at 450. The e-folding time Ψ (the time for the initial post-Flood accumulation rate to drop to 37% of its initial value) was set to 255 years. This produced the accumulation curve shown in Figure 6. The ice softness A was set to 1.0×10^-16 Pa^-3 year^-1 (Bueler, 2016, p. 16). The density of the ice was taken to be 910 kg/m³, and the acceleration due to gravity was 9.81 m/s². Note that, although we are using parameters appropriate to Dome C, this exercise is not a realistic simulation of the real Dome C. Real ice sheets are not isothermal (Paterson, 1980, p. 39), and we have not taken into account possible basal sliding or isostatic sinking (Parrenin et al., 2007).

            The simulation was run with a starting dome radius R₀ of 1000 kilometers, with 101 grid points in both the x and y directions. Because the horizontal plotting surface was 4000 km by 4000 km, this resulted in a horizontal grid resolution of Δx = Δy = 40 km. The large radius was chosen because, in Oard’s Ice Age model, snow is expected to be deposited over large areas in ‘snowblitzes’ (Oard, 2006, pp. 77-81). Thus one would expect accumulation to cover large areas even after just one year of deposition. For each iteration, the ice was allowed to fall for all values of x and y such that x² + y² was less than or equal to the fraction f = 0.85 of the dome radius for the current iteration. There were two reasons for doing this. First, one would expect continuing accumulation to occur everywhere, except perhaps at the periphery of the ice sheet. Second, and more important, the discrete nature of the horizontal grid implies that the number of pixels along the dome radius is slightly greater if the radius is parallel to a horizontal axis than at an angle to it. Because the new value of the dome radius was measured along the x or y axis at each time step, this pixilation effect would cause ice deposited along the axes to fall outside the current dome radius if f were set to 1.0. This would result in the dome radius becoming larger and larger, simply due to roundoff error.

            The results are shown in Figure 7. Initially the ice sheet is a perfectly flat circular slab with a starting height equal to the first year’s total ice accumulation of 12.8084 meters. Because the slab surface is flat (surface slope of zero) the driving stresses are zero, except at the edges. These stresses cause the slab edge to become rounded as the ice spreads radially outward. After 500 years, the ice divide is 2821 meters high, but much of the slab is still flat. Hence, forming a thick Antarctic ice sheet in the time after the Flood seems feasible, provided that the starting post-Flood accumulation is on the order of 10 meters/year, as suggested by Oard (2021). After 1000 years, the dome has increased in height to 3229 meters and become a little broader, but the top is still truncated. At 2500 years, the dome is 3336 meters high, with a slightly truncated top. After  4500 years, the height has subsided a little to 3277 meters, but the dome has become slightly broader with a pointed top (Figure 8).

Consistency Checks

My IDL code did not assume the ice dome to be symmetrical about the lines x = 0 and y = 0, yet one would expect this since precipitation was deposited in circular sheets centered on the origin. That this is the case is encouraging. Also, because the ice is incompressible, the total volume of ice deposited in 4500 years should equal the final dome volume. Summing the total volume of ice deposited over 4500 years and directly calculating the dome’s final volume both yield 8.590×10¹⁵ m³ of ice, as expected (however, when doing this calculation, one must take into account the discrete nature of the ice “blocks”). Finally, it is encouraging that the final dome has a pointed “peak”, with an apparently non-zero slope (Figure 8), as this agrees with  field observations (Vialov 1958, cited in Paterson 1980, p. 22). Note however, that the if one were to closely “zoom into” the divide region, the surface slope would indeed be zero at the divide location, although this is in not apparent in the large-scale depiction in Figure 8. In fact, a zero slope at x = 0 is demanded by symmetry: since the horizontal movement of ice away from the divide depends largely upon the surface slope α, and since u must be zero for x = 0 (by symmetry), the surface slope α must also be zero at x = 0. However, the extent of this “flat” region surrounding the ice divide is too small (compared to the horizontal pixel size) to be seen in my figures.

Discussion

The results are highly sensitive to the choices for Φ and Ψ. If the choices of Φ and Ψ are too low, the final dome height is relatively small, considerably less than 3000 meters. If the choices are too large, the dome becomes extremely thick, with a final height of more than 4,000 meters. After much trial and error, I was able to obtain results similar to Vardiman’s by using the stated choices of Φ and Ψ. Figure 9 is a comparison of the results from the two models. The actual height of the Dome C ice core is 3273 meters (Parrenin 2007, p. 244). After 4500 years, Vardiman’s model yields a final height of 3266 meters, and Mahaffy’s model, combined with Vardiman’s accumulation model, yields a final height of 3277 meters.  

Caveats and Future Research

Biblical creationists should be especially interested in ice behavior at the center of the ice sheet, as ice cores are drilled near the center. Specifically, we would like to know the true thicknesses of annual layers at the center of the ice sheet. Unfortunately, Mahaffy’s model ignores normal stresses/strains, but these cannot be ignored at the center of the ice sheet or at the edges (Paterson, 1980, p. 43). In fact, because it is the surface slope that drives the movement of the ice, the ice at the center of the sheet does not thin at all until the slope at the divide becomes non-zero. Since layer thicknesses depend on normal stresses/strains at the divide location, it may be problematic to use Mahaffey’s model to find these thicknesses. However, this may not be a large source of error, as Weertman (1961, p. 953) noted that these normal or longitudinal strains are only important for smaller ice sheets, about 30 kilometers in width. Since the Antarctic and Greenland ice sheets are much larger than this, it may still be possible to use Mahaffey’s model to obtain estimates of annual layer thicknesses at the center of a large ice sheet. Mahaffey’s model may be modified to take into account variations in topography of the underlying bedrock, and, to a lesser extent, may take into account variations in ice temperature, by “weighting” the value of A(T) toward temperatures near the base of the ice, since this is where the greatest amount of horizontal flow occurs (Figure 1).

The next logical step is to solve the simple cases of isothermal (infinitely) long ice ridges and isothermal symmetrical ice domes, using the more rigorous Blatter-Pattyn (Blatter 1995; Pattyn 2003) solution, which does not ignore longitudinal stresses. Those results can then be compared with those presented here. A detailed discussion of the Blatter-Pattyn model and its solution is provided at http://websrv.cs.umt.edu/isis/index.php/Blatter-Pattyn_model. Because of symmetry, these simplified problems are two-dimensional, which should help reduce computation times. The freely-available Community Ice Sheet Model (or CISM, Lipscomb et al. 2019) uses the Blatter-Pattyn method, but it is a complicated, modular code, and I have had difficulty getting it to work. Hence, it may be necessary to write my own code to obtain these solutions.

 

 

 

Figure Captions

Figure 1. An idealized ice divide. Ice left of the divide flows to the left, and ice right of the divide flows to the right. Ice at the divide location must move straight down. Image courtesy of Michael J. Oard.

Figure 2. A small cubical parcel of ice.

Figure 3. The difference in normal stresses on two opposite cube faces tends to compress or extend the cube, and differences in shear stresses between opposite faces tend to both deform the cube.

Figure 4. Coordinate system used in Mahaffy’s ice sheet model.

Figure 5. Horizontal fluxes (x-direction only) of ice flowing into and out of a column of ice.

Figure 6. The accumulation rate (meters of ice per year) used in this exercise, calculated from Vardiman’s model using values of Φ = 450, Ψ = 255 years, and λ_H/τ = 0.0284 meters of ice per year. Although hard to see on this graph, for large times, the accumulation does not go to zero but asymptotically approaches λ_H/τ.

Figure 7. The ice dome calculated from Mahaffy’s model using Bueler’s diffusive solution for a) 500 years, b) 1000 years, c) 2500 years, and d) 4500 years. The dome height is greatly exaggerated compared to the dome radius.

Figure 8. Side view of ice dome at t = 4500 years.

Figure 9. Height of the ice divide H as a function of time, calculated from both Mahaffy’s and Vardiman’s models and using in both cases the values of Φ, Ψ, and λ_H/τ specified in the text and Figure 6.

References

Alley, R.B. 1988. Concerning the deposition and diagenesis of strata in polar firn. Journal of Glaciology 34(118): 283-290.

 

Alley, R.B., and B.R. Koci. 1988. Ice-core analysis at Site A, Greenland: Preliminary results. Annals of Glaciology 10: 1-4.

 

Alley, R.B., C.A. Shuman, D.A. Meese, A.J. Gow, K.C. Taylor, K.M. Cuffey, J.J. Fitzpatrick, P.M. Grootes, G.A. Zielinski, M. Ram, G. Spinelli, and B. Elder. 1997. Visual-stratigraphic dating of the GISP2 ice core: Basics, reproducibility, and application. Journal of Geophysical Research 102 (C12): 26,367-26,381.

 

Anonymous. Dating a core. Australian Antarctic Division, Department of the Environment and Energy. Posted on antarctica.gov.au. Revised July 25, 2020. (accessed February 23, 2021). [https://www.antarctica.gov.au/about-antarctica/weather-and-climate/climate-change/ice-cores/dating-a-core/]

 

Anonymous. Dating by forward and inverse modelling. Centre for Ice and Climate: Niels Bohr Institute. Posted on iceandclimate.nbi.ku.dk. (accessed February 23, 2021). [https://www.iceandclimate.nbi.ku.dk/research/strat_dating/forward_inv_modelling/]

 

Blatter, H. 1995. Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients. Journal of Glaciology 41 (138): 333 – 344.

 

Bueler, E. 2016. Numerical modelling of glaciers, ice sheets, and ice shelves. Notes for the International Summer School in Glaciology. McCarthy, Alaska. [https://glaciers.gi.alaska.edu/sites/default/files/notes-bueler-2016.pdf]

 

Cuffey, K.M., and W.S.B. Paterson. 2010. The Physics of Glaciers (4^th ed.). Butterworth-Heinemann (an imprint of Elsevier), Amsterdam, Netherlands.

 

Dansgaard, W., and S.J. Johnsen. 1969. A flow model and a time scale for the ice core from Camp Century, Greenland. Journal of Glaciology 8(53): 215-223.

 

Glen, J. W. 1955. The creep of polycrystalline ice. Proceedings of the Royal Society of London, Series A: 519-538.

 

Hammer, C. U., H. B. Clausen, W. Dansgaard, N. Gundestrup, S. J. Johnsen, and N. Reeh. 1978. Dating of Greenland ice cores by flow models, isotopes, volcanic debris, and continental dust. Journal of Glaciology 20 (82): 3-26.

 

Hebert, J. 2014. Ice cores, seafloor sediments, and the age of the Earth, Part 2. Acts & Facts 43(7).

 

Hebert, J. 2018a. Bill Nye, PBS highlight young-Earth evidence. Creation Science Update. Posted on ICR.org. April 27, 2018. (accessed July 7, 2020). https://www.icr.org/article/bill-nye-PBS-young-earth

 

Hebert, J. 2021. Using Vardiman’s young-Earth ice sheet model and a simple computer code to estimate annual layer thicknesses. Creation Research Society Quarterly 57: 175-185.

 

Lipscomb, W. H., S. F. Price, M. J. Hoffman, G. R. Leguy, A. R. Bennett, S. L. Bradley, K. J. Evans, J. G. Fyke, J. H. Kennedy, M. Perego, D. M. Ranken, W. J. Sacks, A. G. Salinger, L. J. Vargo, and P. H. Worley. 2019. Description and evaluation of the Community Ice Sheet Model (CISM) v2.1. Geoscientific Model Development 12: 387-424.

 

Mahaffy, M. W. 1976. A three-dimensional numerical model of ice sheets: Tests on the Barnes ice cap, Northwest Territories. Journal of Geophysical Research 81 (6): 1059-1066.

 

Meese, D. A., A. J. Gow, R. B. Alley, G. A. Zielinski, P. M. Grootes, M. Ram, K. C. Taylor, P. A. Mayewski, and J. F. Bolzan. 1997. The Greenland Ice Sheet Project 2 depth-age scale: Methods and results. Journal of Geophysical Research 102 (C12): 26411-26423.

 

Nye, J. F. 1952. The mechanics of glacier flow. Journal of Glaciology 2 (11): 103-107.

 

Nye, J. F. 1957. The distribution of stress and velocity in glaciers and ice sheets. Proceedings of the Royal Society of London, Series A: 477-489.

 

Oard, M.J. 2005. The Frozen Record: Examining the Ice Core History of the Greenland and Antarctic Ice Sheets. Institute for Creation Research, Santee, CA.

 

Oard, M. J. 2006. Frozen in Time: Woolly Mammoths, the Ice Age, and the Biblical Key to Their Secrets. (2^nd printing). Master Books, Green Forest, AR.

 

Oard, M. J. 2021. Ice core oscillations and abrupt climate changes: part 3 – Large-scale oscillations within biblical history. Journal of Creation 35 (2): 107-115.

 

Parrenin, F., G. Dreyfus, G. Durand, S. Fujita, O. Gagliardini, F. Gillet, J. Jouzel, K. Kawamura, N. Lhomme, V. Masson-Delmotte, C. Ritz, J. Schwander, H. Shoji, R. Uemura, O. Watanabe, and N. Yoshida. 2007. 1-D-ice flow modelling at EPICA Dome C and Dome Fuji, East Antarctica. Climates of the Past 3: 243-259.

 

Paterson, W.S.B. 1980. Ice sheets and ice shelves. In: Colbeck, Samuel C. (editor). 2012. Dynamics of Snow and Ice Masses. Academic Press, New York, NY.

 

Pattyn, F. 2003. A new three-dimensional higher-order thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes. Journal of Geophysical Research 108, No. B8. doi:10.1029/2002JB002329.

 

Vardiman, L. 1993. Ice Cores and the Age of the Earth. Institute for Creation Research. San Diego, CA.

 

Vardiman, L. 1994. An analytic young-Earth flow model of ice sheet formation during the “Ice Age”. In: Proceedings of the Third International Conference on Creationism. Walsh, Robert  (editor). Creation Science Fellowship, Pittsburgh, PA, pp. 561-568.

 

Vardiman, L. 2001. Climates Before and After the Genesis Flood: Numerical Models and Their Implications. Institute for Creation Research, El Cajon, CA, pp. 41-68.

 

Vialov, S. S. 1958. Regularities of glacial shields movement and the theory of plastic viscous flow. IASH Publication 47 (Symposium at Chamonix, 1958  - Physics of the Movement of Ice): 266-275.

 

Weertman, J. 1961. Equilibrium profile of ice caps. Journal of Glaciology 3 (30): 953-964.