We have undertaken a simplified calculation of orbital forcing back through the Cretaceous to the Late to Middle Jurassic from 65 to 190 Ma. So long as the Earth has a continental ice volume, orbital forcing will impose a 400-ky periodicity upon glacioeustasy and thereby on fourth-order sequence stratigraphy cycles. Similarly, third-order cycles were defined by orbital forcing of 2.4 ± 0.4 my (predominately 2.0- and 2.8-my cycles). These concepts greatly simplified the task of unraveling sequence stratigraphy. Our sea-level calculations are comparable with stratigraphic observations and the results are consistent with a persistent continental ice volume throughout the Late to Middle Jurassic and Cretaceous. In general, they compare well with the Arabian Plate Maximum Flooding Surfaces and the Cretaceous and Jurassic stage boundaries, within the limits of the recognized stratigraphic time scales. We used simple Parametric Forward Models (PFMs) for modeling changes in sea level, subsidence, and sedimentation and noted that PFMs can be applied to other tasks. The results will provide for rapid, cost-effective forward modeling on tasks such as reservoir characterization and flow simulation.


Sea level changes with time, both daily—diurnal tides—and over geologic time—principally through cyclic climatic changes. Changes in sea level affects marine sedimentation and hence stratigraphy. Orbital forcing is the imposition of cyclicity, due to variations in the Earth’s orbit, on the stratigraphic record. The mechanism behind this cyclicity is climatic change driven by well-understood and predictable changes in the parameters of the Earth’s orbit around the Sun. We consider the following sequence of events:

  • Orbital forcing controls climate (temperature).

  • This has a direct effect on the making and melting of continental ice.

  • This results in changes to the volume of the oceans.

  • This causes variations in sea level and in the types of sediment deposited.

So long as the Earth has an ice volume upon a continent, periodicities in the Earth’s orbit will produce periodicities in sea-level fluctuation. In favorable circumstances, the climatic changes are imprinted on the stratigraphic record in cycles of about 20,000 thousand years (20 ky) to 400 ky years duration. In addition, modulation occurs when two or more sine waves beat together. Modulation introduces important periodicities that have durations of several million years.

Because the readership of this paper might range from astrophysicists to geologists, we have written it in four parts. Part 1 provides an introduction to Orbital Forcing. Part 2 presents the calculation of Orbital Forcing in a form relevant to this paper. Part 3 deals with the effect of Orbital Forcing on the Earth’s climate and sea level. Part 4 considers the resultant relationship between Orbital Forcing and Sequence Stratigraphy.

We acknowledge that numerous processes affect climate, the volume of the ocean basins, and the volume of water in them. Figure 1 shows the effects of periodic and chaotic events on the global ocean volume and climate. Chaotic events are those that occur randomly, or are of very long period, or are poorly understood. In contrast, orbital forcing is periodic and will produce a cyclicity in the climate and hence a cyclic stratigraphic record. However, chaotic and periodic processes may interact. For example, note that orbital forcing has to be combined with the favorable disposition of continents at or around the Poles (as at the present time) for widespread continental glaciation to take place, and hence for sea level to be affected. The corollary is that so long as the Earth has an ice volume upon a continent, periodicities in the Earth’s orbit will produce periodicities in sea-level fluctuation that can be recognized in the stratigraphic record.

Sequence stratigraphy has become an extremely important tool in petroleum exploration and production but there is always room for improvement. We seek to bring orbital forcing of glacioeustasy (Hays et al., 1976; Zachos et al., 2001) more prominently into sequence stratigraphy thinking by examining its role in third- and higher-order sequence stratigraphic cycles from the Cretaceous to mid-Jurassic (from 65 Ma to 190 Ma). Calculations of orbital forcing are used to compute sea level as a function of Antarctic continental ice volume and the results are compared with key observations of the sequence stratigraphy of the Arabian Plate by Sharland et al. (2001) and of Europe (de Graciansky et al., 1998).

When geologists try to apply the results of astronomers to the question of orbital forcing on the scale of a few hundred million years, there is a paradox. On the one hand, the astronomers are confident that the major periodicities observed in the last few million years should prevail throughout the Phanerozoic (Berger et al., 1992). On the other hand, they note an increasing uncertainty among orbital components as we move to time scales beyond 50 Ma (Laskar, 1990, 1999). Part of the problem arises from orbital resonances that appear and disappear in a chaotic manner. Furthermore, even small errors in starting conditions render formal calculations to be substantially in error, from the astronomer’s perspective, on the time scales of interest here (Laskar, 1995, p. 20). We propose a middle ground that may be less than satisfactory to astronomers, but useful to problem solving in sequence stratigraphy. To this end, we first put forward a simple overview of orbital forcing. For the convenience of readers who are remote from Orbital Forcing and/or Parametric Forward Modeling, technical words are shown in bold face when defined in the nearby text.


Orbital Variations

The hypothesis that climatic change was linked to periodic changes in the Earth’s orbit due to the influences of other bodies in the Solar System (Figure 2) was first proposed by Croll in the mid-nineteenth century. It was refined by Milankovitch (for more information see Berger and Loutre, 1991; Hays et al., 1976; Murray and Dermott, 1999) and his name is now linked to the astronomical theory of cyclic climate variations.

These periodic (cyclic) perturbations in the Earth’s orbit (Table 1) are the results of the gravitational effects of planetary bodies on the Earth. Their combined influence causes continuous variations in the seasonal and geographic distribution of solar energy (or insolation) received at the Earth’s surface and thus acts as an external forcing mechanism on the Earth’s climate.

Table 2 and Figure 3 provide an Earth’s eye view of the Solar System. A quick review of this data provides insight into the major parameters of the Solar System. Note that the largest planet (Jupiter is about 300 times the mass of Earth, but it is only approximately 1/1000 the mass of the Sun. Venus and Earth are of comparable mass; Mars has about 1/10 the mass of the Earth. The relative gravitational effect of each planet on Earth is proportional to the mass of the planet and inversely proportional to the square of the distance between the planet and Earth. Jupiter has the largest effect on the Earth’s orbit; the effect of Venus is about half that of Jupiter and the effect of Saturn, about 1/10 that of Venus; and so on. Note that there is a large jump in distance from the Sun beyond Saturn (planet six) and that for purposes of calculations relevant to Earth history on the scale of hundreds of millions of years, the effects of Uranus, Neptune, and Pluto can be ignored.

In the following commentary, we will be interested in the period of various terms taken alone and in the modulation created by the interaction of two or more terms. Eccentricity, tilt, and precession can all be written



where E is eccentricity, tilt, or precession, A is an amplitude for each term, (1/p) is the frequency (the inverse of period) of each term, t is time from starting conditions (in our case, back in time), and phs is a phase angle for starting conditions.

As an example of simple modulation, consider two prominent terms of eccentricity, 95 ky and 99 ky with amplitudes of 0.8 and 0.6, respectively. Exactly what they measure is not important at this point. Consider both to be going around the Sun and the faster 95 ky term to be ‘chasing’ and periodically ‘overtaking’ the slower 99-ky term. The time required for 95 ky to lap 99 ky is given by subtraction of the two frequencies. Converting back to period, we find the modulation to have a period of 2.34577 million years (2.35 my, and within rounding error of the actual modulation of eccentricity, 2.4 my). Where the two sine waves cancel each other, the modulation forms a node. Where the two sine waves reinforce each other, the modulation forms an antinode. The amplitude at nodes is given by subtraction (0.8 - 0.6 = 0.2); the amplitude at antinodes, by addition (0.8 + 0.6 = 1.4). In the real world, modulation involves more than two periods—a point to which we will return.


The eccentricity of each planet’s orbit is a measure of its ellipticity. It is measured as the difference between perihelion (the point in the planet’s orbit closest to the Sun) and aphelion (the point in the planet’s orbit farthest from the Sun) divided by the average of the two distances (Figure 4a). The orbits of the various planets vary from nearly circular to strongly elliptical. Venus and Earth have very low eccentricity whereas Mercury has an eccentricity of the order of 40 percent. Mars, Jupiter, and Saturn have eccentricities of about 10 percent at present. The eccentricity of the Earth’s orbit changes with time from almost zero to about six percent—the longest term is 400 ky and there are numerous terms of about 100 ky. Important modulation occurs at intervals of 400 ky and 2.4 million years (my).

Given the length of the year of the various planets, it is easy to see how the orbits become eccentric. Consider, for example, the gravitational interactions between Venus, Earth, and Jupiter as they orbit the Sun. Jupiter takes about 12 Earth years to make one revolution about the Sun; Earth, one year; and Venus, 0.62 years. These differences in length of year dictate that Earth has the opportunity to see Jupiter, Venus, and the Sun in the following three very different configurations.

  • When Earth crosses the line connecting Jupiter and the Sun, and Venus is on the other side of the Sun, the effect of Jupiter is to pull Earth farther away from the Sun.

  • Conversely, when Earth crosses the line connecting Venus and the Sun, and Jupiter is on the other side of the Sun, the effect of Venus is to pull the Earth closer to the Sun.

  • When Jupiter, Earth, and Venus are all on the same line to the Sun, the tug-of-war between Jupiter and Venus goes somewhat in favor of Jupiter.

Thus, we can begin to see the problem facing the planetary astronomers and the vast difference between their time scales and those of the geologists. Earth and Jupiter interact in some fashion each year and the location of Venus at that time is a separate question. The problem generalizes to all planets having some effect on all other planets during some part of their annual orbital cycle. The net result is an Earth that has an elliptical orbit that is also eccentric, as the values of perihelion and aphelion change in a cyclic manner. The main parameter in orbital forcing is eccentricity, with tilt and precession having only a minor part to play.

Tilt and Precession

Tilt and Precession are variations in the angle and direction of the Earth’s axis.

Variations in the Obliquity or Tilt of the Earth’s rotational axis relative to its orbital plane result in stronger or weaker differentiation of the seasons—if the Earth’s axis were perpendicular to the orbital plane, there would be virtually no seasons. The tilt angle varies between 21.1° and 24.5° and is presently 23.4° (Figure 4b). Obliquity includes terms with periods of about 40 ky and 54 ky and shows complex, low amplitude 1.0 to 1.4 my modulation.

Precession is the gradual change in the direction in which the Earth’s axis is pointing. Like a spinning toy top, its axis of rotation executes a slow precession (Figure 4c). It is caused by gravitational torque exerted on the spinning Earth by the Moon and Sun. The main effect of precession of the Earth/Moon System is in combination with eccentricity. The warmest summers occur when summer solstice is at perihelion during times of high eccentricity (Figure 5).

In our calculations of sea level (below), the frequencies of the classic periodicities of precession and tilt (19,000–23,000 years and 41,000–54,000 years) are too short to produce large changes in sea level during the Mesozoic. Nevertheless, we deal with them for completeness of our calculations. Precession of the Earth/Moon System basically samples the much longer-term modulation of eccentricity (producing Climatic Precession), and tilt has a complex modulation of its own on the order of 1.0 to 1.4 my, with alternative larger and smaller antinodes. Additionally, the periodicities of both precession and tilt change systematically back through time with changes in the Earth/Moon precessional constant, k. Precession and tilt are fully taken into account in our calculation of insolation values that drive our sea-level calculations.

For the last few hundreds of thousands of years, the effect of precession and tilt appears to center around summer melting of Northern Hemisphere ice sheets when precession and/or tilt maxima has coincided with the Northern Hemisphere summer solstice at perihelion. Stated more generally for longer time scales, as indicated in Figure 5, we propose extended periods of relatively low eccentricity are times when the Antarctic ice sheet continues to get larger and larger. Times of relatively high eccentricity, combined with precession and tilt maxima at times of Southern Hemisphere summer solstice at perihelion tend to melt ice over time scales on the order of a few hundred thousand years; commonly the warming half of 400-ky cycles.


In order to examine the effects of orbital forcing on sequence stratigraphy we have calculated a simple version of orbital forcing from 65 Ma back to 190 Ma and examined the origins of periodicity in both amplitude and amplitude modulation. Figure 6 summarizes our overall research strategy. We first calculated an approximation of the variations in the Earth’s orbit back through time (Figure 6a). In Part 3, we have calculated orbital forcing of sea-level history at 1,000-year time steps (Figure 6b) with a parameterization such that our calculated sea-level file converges with stratigraphic observations (Figure 6c). Thus calibrated, Part 4 compares our sea-level file with other stratigraphic observations.

Mechanics of Orbital Forcing

The high-precision approach to the mechanics of orbital forcing is through the numerical integration of the Solar System. Starting with the present configuration of the planets around the Sun, the system is stepped back in time in small increments. The gravitational attraction of each mass on all other masses is calculated at each step to determine the position of each planet at the next step back in time (Quinn et al., 1991; Laskar, 1990). The high-precision calculation is extremely computer-intensive, time-steps being on the order of a few Earth-days. As a result of Fourier analysis, the results can be encapsulated in the fundamental frequencies,g and s, for each planet (columns to the right in Table 2).

The g value for each planet describes the time required for the semi-major axis of the elliptical orbit to travel one orbit or 360° around the Sun. It is expressed as a frequency (rate). The fundamental frequency, s, is a similar representation for the cyclical relationship of each planet’s orbit to the plane of mass of the Solar System. For practical purposes, the plane of mass of the Solar System is the orbit of Jupiter about the Sun. The orbit of each planet periodically makes a slight angle with the plane of mass. The frequency of this up and down motion is the fundamental frequency, s. Much of the variation in tilt is due to s.

Relying on g and s values for each planet, an approximation of Orbital Forcing can be made from simple trigonometric series as expressed in Equation 1. This is computationally simple as time-steps (here, 1 ky) are a matter of sampling known sine waves, and the calculations can be run on a high-end personal computer. A further advantage of this approach is the ease in which the role of each term, or combination of terms, can be examined in isolation. The complex becomes simple when taken apart.

Our eccentricity calculation was based on the Matthews et al. (1997) 1,000-year time-step parametric approximation (hereafter, MFD97). When MFD97 is plotted against the more rigorous Laskar (1990) time series calculation (LA90), the relationship is slightly non-linear. We applied a linear correction to both ends of the MFD97 file to make the relation to LA90 linear and to bring minimum and maximum values within bounds of the LA90 target data. The resulting file (this paper) has a correlation coefficient of 0.95 with regards peak for peak amplitude comparison to LA90 target data.

Simplified estimates of tilt and climatic precession were calculated in a manner similar to the MFD97 eccentricity calculations. The value of k (precession of the Earth/Moon system) was increased by 0.066 arc seconds/year every 2.4 my, using an initial k value of 50.4712 (Berger et al., 1992). The resulting estimate of tilt angle has a correlation coefficient of 0.94 with regard to peak-for-peak comparison to LA90 target data. The resulting estimate of precession has a correlation coefficient of 0.984 with regards peak-for-peak amplitude comparison to LA90 target data. From climatic precession, the derivation of the angle omega (the longitude of perihelion; expressed as an angle relative to the Northern Hemisphere spring equinox) is straightforward. In all of our calculations, age discrepancies between our calculations (this paper) and LA90 are minor for the present exercise.

From our files of eccentricity, tilt, and omega, the solar insolation values for any month at any latitude can be calculated in time steps of 1,000 years or greater. We utilized insolation files for July 30°N and for January 70°S to drive a simple PFM of glacioeustatic sea-level fluctuations.

Table 3 presents rank, components, amplitude, phase, frequency, and period for the 19 terms of the simplified eccentricity calculations of Matthews et al. (1997). The components in Table 3 are the combinations of fundamental frequencies relevant to planets 1 to 6, ranked by amplitude derived from comparison with the high-precision calculation of eccentricity (LA90). Consistent with our qualitative discussion of Table 2, the largest term in defining the eccentricity of the Earth is its interaction with Jupiter and Venus and the resultant 400-ky cycle is widespread in our sea-level calculations. It is significant that g5 (Jupiter) or g2 (Venus) are involved in terms 2 to 5 and these are the four strongest ‘100-ky’ terms. Also of importance is term 8, (g4-g3) that has a period of about 2.4-my. Although this term is only about one-quarter of the amplitude of term 1, combinations produce strong (g4-g3) modulation, a point to which we will return.

Orbital Cycles: a ‘Tuning Fork’ for Geologists

Our analysis of orbital forcing begins with a 400-ky pure sine wave (Figure 7a). Note that this is the first term in the eccentricity trigonometric series presented in Table 3. This is perhaps the closest to a ‘tuning fork’ for geologic time that geologists will ever find.

Figures 7b and 7c demonstrate 400-ky modulation in combination with Mars (terms 2 and 3; 95 and 124-ky) and Earth (terms 4 and 5; 99 and 131 ky). When considering the combined effect of two terms, such as 2 and 3 (Table 3), the frequency of modulation is given by the subtraction of term frequencies, here (g4-g5) - (g4-g2) = (g2-g5), which is the 400-ky periodicity. The amplitude of this modulation at node is given by the subtraction of amplitude 3 from amplitude 2; and at antinode, by the sum of the two amplitudes. Note that nodes in modulation in Figures 7b and 7c correspond to troughs on the sine wave Figure 7a and that antinodes in Figures 7b and 7c correspond to peaks on the sine wave in Figure 7a. Figures 7d and 7e go on to relate these nodes and antinodes to the sum of all ‘100-ky’ terms in Table 3 and to eccentricity. Modulation at 400-ky is clearly represented in both Figures 7d and 7e, a point to which we will return.

Figure 8 presents a similar analysis with regard to (g4-g3), the 2.4-my cycle. As with Figure 7a, it begins with a pure sine wave that persists throughout at least the Phanerozoic as a ‘tuning fork’ of geologic time. However, it will be seen that the amplitude of term 8 in Table 3 is only about 25 percent of that of term 1 (Table 3 and Figure 7).

Figures 7b and 7c and 8b and 8c involve modulation of the four ‘100-ky’ terms (2, 3, 4, and 5 of Table 2), but in different combinations. It appears that the modulation is much more spectacular in Figures 8b and 8c than in Figures 7b and 7c. This is because g3 and g4 (Table 2) are closer together than any other g values. This produces a longer period modulation of the ‘100-ky’ signal, which is visually more appealing than the shorter period modulation in Figures 7b and 7c. In fact, the node/antinode ratio for Figures 7b and 7c is 12 percent whereas for Figures 8b and 8c it is 15.5 percent, a point to which we will return.

With regard to Figure 8, nodes and antinodes show up well in the sum of all ‘100-ky’ terms (8d) and in eccentricity (8e). Indeed, from examination of the modulations depicted in Figures 7 and 8 for all combinations of terms in Table 3 that reduce to the requisite (g2-g5) [400-ky] and (g4-g3) [2.4-my] modulation, the same relation of modulation to term 1 or term 8 sine waves prevails. This is a phase lock in as much as anything that affects the phase of terms 1 or 8 will have a similar affect on the phase of the related modulation. Note that this result supplants the related approximate calculations of Matthews et al. (1997).

We examined further this important assertion in two steps. First, we inspected the MFD97 eccentricity calculation for phase lock of the 2.4-my terms (Figure 8) back to 190 Ma. Not surprisingly, the phase lock remained in force. Second, we undertook alternative calculations of eccentricity in which g and s values change in accordance with Laskar (1990, Figures 8 and 9). His results suggest small, systematic variations in g and s values over the last 200 my.

When we introduced these changes in g and s and compared the results to MFD97 (with its constant values for g and s), we observed phase shifts in the sum of ‘100-ky’ terms, but the phase lock prevailed. Further, we counted the antinodes from the present back into the Jurassic and correlated them from one calculation to another. We found that recognizable ‘signatures’ occur within the 2.4-my modulation of the sum of ‘100-ky’ terms. In Figure 8d, for example, the antinode at 3 to 4 Ma has a different appearance from the one at 5 to 7 Ma. Even with sizeable phase shifts resulting from slight changes in g and s values, the ‘signature’ of each particular antinode (counting back in time from the present) commonly remained intact. This result is consistent with the calculations of Matthews et al. (1997).

Thus, we see that the regular, cyclic behavior of orbital parameters (eccentricity, tilt, and omega) produces changes in solar insolation that will ultimately affect continental ice volume and therefore the ocean water volume and sea level, and finally the stratigraphy. These changes include cycles with periods of about 400 ky and 2.4 my. These, and related phase-locked modulations, will persist even if there are chaotic shifts in the phases relating the various orbital parameters. Therefore, we can expect orbital forcing of sequence stratigraphy even over time scales where the astronomers suspect phase shifts are certain to have occurred.


Changes in the cryosphere (the part of the Earth that is covered by permanent ice) directly affect the volume of water in the oceans and hence global sea level. The changes are effected by the periodic (cyclic) perturbations in the Earth’s orbit that cause variations in solar insolation and hence the climate. As an example, consider the Quaternary. For the last 2 my or so, high-latitude winters have been cold enough to allow snow to accumulate, but it was only when summers were cold that the snows of previous winters did not melt completely. This occurred when the Earth was in aphelion and favorable phases of orbital eccentricity, tilt, and precession combined. When this process continued for long enough, ice sheets began to form. The process was reinforced by positive feedback effects due to the snow and ice reflecting increasing amounts of solar radiation, and by changes in the distribution of ocean currents that convey heat from one part of the Earth to another. Conversely, major interglacial periods occurred when eccentricity, tilt, and precession gave the Northern Hemisphere the greatest amount of summer insolation.

It is unlikely that Northern Hemisphere glaciers of any significant size occurred during the Mesozoic. Thus attention shifts to the Southern Hemisphere. Whereas during the Pleistocene ice age the Northern Hemisphere had an ocean at the pole and several continents at various distances from the pole, Antarctica presented a much simpler continental-centered geometry throughout the Cretaceous and much of the Jurassic. We constructed a simple Parametric Forward Model (PFM) to calculate the size of an Antarctic continental ice sheet, and hence the related change in water volume in the ocean and thereby sea level, back to 190 Ma. We then related the results to Mesozoic sequence stratigraphy (Part 4).

A Simple Parametric Forward Model (PFM) for Mesozoic Sea Level

Basic Concepts of Simple Parametric Forward Modeling

In forward modeling, reliance on the understanding of geological processes allows high-resolution prediction of various geological parameters. We employed our forward model to relate sea level to solar insolation values derived from the orbital calculations described above. Inasmuch as the PFM approach is quite different from process modeling, we briefly review the basics and then go on to describe our sea-level PFM and resultant calculations.

In this paper, we are attempting to model the outcome (glacioeustatic sea-level fluctuations) of Atmosphere/Hydrosphere/Cryosphere interactions from 65 Ma back to 190 Ma. To approach this problem in strict dynamic process modeling mode would be an enormous task. Such models for atmosphere and ocean circulation typically involve the transfer of physical properties among thousands of cells by solving simultaneous nonlinear differential equations at small time steps. The models are mostly run in perpetual mode (seasonal or annual) that produces equilibrium snapshots of output under variously specified boundary conditions—for example, the work of CLIMAP (1976, 1981). Time series dynamics is commonly reconstructed by comparing a few of these snapshots (see COHMAP, 1988; Webb et al., 1993). When these models are run in time-series mode, the goal is generally on the scale of the dynamics of the annual cycle (Oglesby, 1982).

As an alternative to dynamic process modeling, we contend that problems can be solved rapidly on distant time scales by employing simple PFMs. One does not need to calculate every single step in the process from ocean to atmosphere to ice sheet to capture the basic idea that warm summers tend to melt more ice. A PFM attempts to replace a full numerical solution of the differential equations describing complex physical, chemical and biological processes, with a step-wise linear time series parametric calculation that can be compared to target data. The basic goal is to keep the PFM as simple as possible and still have it converge on the target data; make the model more complex only if initial comparisons are unacceptable. For a PFM to be acceptable there must exist at least one set of model parameters such that the model and the target data are sufficiently in agreement and the parameters are geologically reasonable.

Note that an acceptable solution to the forward problem is not a unique solution but demonstrates the capability of the model to converge with the data. There remains the inverse problem (Matthews and Frohlich, 1998, p. 377); namely, the finding of all sets of model parameters such that (1) the model and the target data are sufficiently in agreement, and (2) the parameters are all geologically reasonable. Indeed, the reader may ask, “How do the authors know that particular parameter?” Commonly, the answer is, “They don’t.” The parameter is simply geologically reasonable. The simple PFM yields an acceptable solution with simultaneous use of that parameter and a number of related parameters. Change one parameter, and the model no longer converges on the target data. Change some other parameter or parameters and the model again converges on the target. The solution to the inverse problem is beyond the scope of this paper. Our first step was to calculate a sea-level file and demonstrate by model/data comparison that it has credibility (the task of this paper).

Model Description

Our PFM sea-level calculation was based on solar insolation for July at 30°N and January at 70°S in 1,000-year time steps. We prescribed a starting sea level (the model rapidly takes over), parameters relating to Antarctic ice-sheet dynamics, and specified maximum and minimum sea levels during each model run. Net ice-make/ice-melt was calculated for each time step as a result of variations in insolation (Figure 9).

Our model for Antarctic ice is simple, yet follows the lead of more detailed process modeling—for example, that of Fastook and Prentice (1994). To begin, there is a specified continental topography capable of accumulating snow year-round. Without such a continental configuration, ice cannot accumulate sufficiently to affect the volume of the oceans. The ice sheet itself is taken to have a prescribed thickness and a variable radius consistent with sea level at each time step. New snow accumulates over the area of specified continental topography or over the calculated area of the ice sheet, whichever is larger. Summer ice melting occurs around the perimeter of the ice sheet, to a prescribed distance in from the ice margin if the 70°S critical insolation parameter is exceeded.

We perceive the physics of January insolation at 70°S to be straightforward. There is a critical value for insolation below which ice continues to accumulate on Antarctica year-round (the situation today), and above which some summer precipitation falls as rain and the ice volume decreases accordingly.

The physics of the Northern Hemisphere insolation is more complex. It is possible that ‘warm saline bottom water’ (more generically ‘reverse thermohaline circulation’) (Brass et al., 1982) may deliver heat to the waters around Antarctica, influencing ice accumulation or ice melting (Prentice and Matthews, 1991). Alternatively, it is possible that orbital forcing effects the deep-sea residence time of normal thermohaline circulation. The longer the average residence time, the more CO2 remains in the deep ocean; the shorter the residence time, the higher the partial pressure of CO2 in the atmosphere, although exact mechanisms and event sequencing may be complex (Boyle, 1988; Broecker and Peng, 1989; Broecker and Henderson, 1998). Regardless of the specific physics, the model proposes that there is a critical value for summer insolation at 30°N below which ice accumulation is increasingly favored and above which ice melting is increasingly favored (Figure 9). Ice-making/ice-melting rates associated with critical insolation values may vary independently or may be tied together as a system.

Sea-level Calculations A and B

Previous computations of sea levels indicated that some boundary conditions generated situations where maximum or minimum ice volumes persisted indefinitely; this produced ‘hang up’. Figures 10 and 11 present sea-level calculations A and B for the period 65 to 190 Ma. The objective of calculation A was to find one set of boundary conditions that does not result in hang up. Calculation B was devised to come closer to realistic conditions for the Jurassic. The orbital-forcing input files were the same for both calculations and differences in the sea-level files were caused by the choice of model parameters.

Calculation A was purposefully set up with boundary conditions that are overly optimistic for the existence of dynamic ice sheets. By specifying an ice-prone continental topography that is probably unrealistically large, ice can once again accumulate even if the glacier becomes relatively small. By setting the critical insolation value at 70°S relatively low, there is generally enough melting to keep the calculation from locking onto maximum ice volume. Finally, allowing sea-level amplitude of 50 m relative to an ice-free world provides a wider range over which to observe the natural tendencies of the calculation. Given these three probable excesses, the calculation does indeed remain within an acceptable range of minus 35 + 15 m (relative to an ice-free world) throughout calculation A.

Calculation B was calibrated to yield an acceptable shape comparison with the Haq et al. (1987) intervals LZA-1 (see Figure 6) and LZA-4.1 through 4.3, and with the Haq et al. time scale adjusted to that of Gradstein et al. (1995). Comparison of the low-resolution Haq et al. coastal onlap curve relative to the detailed sea-level prediction renders this calibration qualitative. These boundary conditions were used to generate the B file from 65 to 190 Ma. The main difference between the A and B calculations (Figures 10 and 11) is that B does hang up around maximum ice-volume conditions from time to time (e.g. 100-105 and 155-157 Ma), but then recovers from these intervals of low sea level.

Comparison of A and B calculations in Figures 10 and 11 provides some indication of the robustness of the ‘third-order signature’. At the scale of these figures, the peaks and valleys result from the 400-ky modulation that is translated into 400-ky periodicity by the asymmetry of ice-making (slow) and ice-melting (fast) in the sea-level calculation. Thus, in looking for a ‘third-order signature’ in the sea-level files, we are looking for the arrangement of 400-ky peaks within a larger time series bounded by major glacioeustatic lowstands.

At least four signatures are noteworthy.

  • A few single peaks that are isolated by at least 2 my on either side stand out clearly in both calculations A and B (e.g. 79, 95, 107, 121, 137.5 and 163.3 Ma).

  • In other places, there are three or more 400-ky peaks arranged in ascending order and then terminating abruptly (e.g. 71–72, 128–29, and 153–155 Ma).

  • Elsewhere, the highest 400-ky peak in a third-order cycle occurs earliest and is followed by lower and lower second, third and fourth peaks (e.g. 76–77, 140.7–142, and 145–147 Ma).

  • Finally, there is in some places an overall rounded shape to the arrangement of four or more 400-ky peaks (e.g. 82–85, 124–127, and 181–183 Ma).

Cycles in Sedimentation

Natural Cycle Order and Amplitude

The Enclosure compares our model results with observational data on sequence stratigraphy. In this section we will consider only the columns on the left in each panel of the Enclosure. Column 1 provides a count of (g2-g5) 400-ky peaks from the present back to 190 Ma. This 400-ky periodicity comes into orbital forcing in two ways. First, it is the highest amplitude pure sine wave in the trigonometric series describing eccentricity (Table 3). Second, it is a prominent modulation of 100-ky eccentricity terms. Asymmetry of processes in the sea-level calculation (ice-making is slow; ice-melting is fast) translates amplitude modulation (such as depicted in Figures 7b, 7c) into periodicity (Figures 6b, 10, 11 and 12). This is the major tuning fork of geologic time. The period of this cycle is precisely defined as (g2-g5), here taken as nominally 400-ky (or 404-ky by current best calculations; Laskar, 1990). We considered this to be the natural definition of ‘fourth-order cycles’.

Similarly, we defined third-order cycles as closely related to the (g4-g3) 2.4-my modulation depicted in Figure 6a, and 8d. We counted antinodes on the (g4-g3) modulation from the present back to 190 Ma and displayed these numbers in column 3 of each panel of the Enclosure.

The numbers in column 3 correspond approximately to whole-number multiples of the fourth-order cycles. The lack of absolute correspondence is because not all third-order cycles are precisely (g4-g3), for at least two reasons. First, as indicated in Figure 8e, there is room within the (g4-g3) node for more than one fourth-order cycle. By our sea-level calculations, maximum glaciation occurs toward the end of those fourth-order cycles that exhibit prolonged relatively minor eccentricity. In Figure 8e, such glacial fourth-order cycles end at 2.4, 4.4, and 7.2 Ma, thus defining 2.0-my and 2.8-my third-order cycles. Second, taking classical geological stage names to be derived from unconformity-bounded shallow-marine sequences, we attempted to point out likely major unconformities by our numbering and color-coding of third-order cycles. In column 3 of each panel of the Enclosure, white indicates unconformities in shallow-water sequences and blue highlights a lowstand immediately preceding a significant rise in sea level. The third-order cycle count is further color coded primarily on the correlation of sea-level calculations A and B with the timing of boundaries, and secondly on the robustness of the third-order ‘signature’ in these two sea-level files. Good agreement is indicated in green, fair agreement in yellow, poor agreement in red.

For those parts of sea-level calculation B that do not involve multiple excursions to minimum glacioeustatic sea level (e.g., 70-85, 115-130, and 145-154 Ma), the typical amplitude of third-order cycles is 17 m and that of fourth-order cycles is 9 m. Where 100-ky (eccentricity) periodicity can be recognized, its amplitude is generally 4 m whereas the ubiquitous 20- to 40-ky (precession and tilt) high-frequency flutter has an amplitude of about 2 m (see Figure 12).

For those intervals of sea-level calculation B that intermittently reach maximum glacioeustatic lowering of sea level (e.g. 93-95, 100-105, and 155-160 Ma), amplitude at various periodicities tends to increase. This is because the basic math of the sea-level calculation involves ice-making and ice-melting volumes that are proportional to the size of the ice sheet. Third-order cycles average about 21 m; fourth-order, about 13 m; 100-ky eccentricity, about 6 m and 20- to 40-ky precession and tilt, about 5 m. Whereas 100-ky peaks and valleys tend to be more recognizable because of higher amplitude, there are still fourth-order cycles over which 100-ky cyclicity is not distinguishable from the background variation in the 20- to 40-ky precession and tilt.

Thus, we can confidently compare third- and fourth-order cyclicity of sea-level calculation B with previous work. However, we caution that for large parts of the Jurassic and Cretaceous, recognition of consistent cyclicity at higher than fourth-order may be an intractable problem that lies not with the stratigraphic sections as recording instruments, but rather with the sea-level history itself. This is because precession and tilt will only be recorded by some attribute of the section that has a faster response time than Mesozoic sea level changes (see for example, the work of Kuhnt et al., 2001).

Nomenclature of Third-order and Fourth-order Cycles

We have employed a decimal nomenclature system to describe the highstand and lowstand of every fourth-order cycle, and this number is displayed in column 5 of each panel of the Enclosure. As an example, consider the Coniacian. The first digit is the third-order cycle count, beginning with the maximum glacioeustatic sea-level lowering preceding the flooding event that begins the third-order cycle. Thus, we named the lowstand at 88.87 Ma as 38.0; the next red dot, 38.1; the next blue dot 38.2; and so forward to the lowstand at 86.85 Ma. This we renamed as 37.0 (rather than 38.10) and resumed counting highstands odd, lowstands even until we reached lowstand 37.10 at 84.83 Ma. Thus, we have given a name (number) to each fourth-order highstand or lowstand. Each blue dot represents a fourth-order cycle boundary and each red dot represents a fourth-order maximum flooding surface (mfs). Where identifiable, we noted third-order maximum flooding surfaces as an ‘MFS’. By our definition, each cycle boundary and flooding surface was the result of glacioeustatic changes.

In columns 6 and 7, we provided an indication of the magnitude of flooding from lowstand to highstand as might be registered in shallow-water, or moderately deep-water, stratigraphic sections. The number indicates the magnitude of sea-level rise. Color-coding indicates the likelihood of this being observed in the geologic record; green, good; yellow, fair; red, poor. The numbers in columns 6 and 7 are estimates of glacioeustatic sea-level rise and should be adjusted for isostatic compensation to arrive at the amount of section that might be represented by these events (see Matthews (1984, p. 134-141, for a simple review.) For relatively deep-water stratigraphic sections, the measurement of the potential new accommodation space is straightforward. We assumed, for example, that lowstand 38.0 in a relatively deep-water stratigraphic section would show a key sedimentary facies that could be recognized as the shallowest water on either side of the interval. Thus, 38.5 would have 19 m (plus isostatic adjustment) more accommodation space than 38.0.

For shallow-water facies, the availability of new accommodation space is further complicated by previous highstands of the sea. For example, we have not indicated any credible flooding event for highstand 38.1. Assume for the moment that shallow-water deposits filled the accommodation space at 39.7. Even though lowstand 38.0 gave way to a considerable sea-level rise, accommodation space would have been filled by deposition at 39.11 or 39.13, leaving little accommodation space to be filled at 38.1. Thus, we did not flag 38.1 as a likely candidate for a fourth-order maximum flooding surface. In contrast, sedimentation to fill accommodation space at highstands 38.7 and 38.9 would still leave 9 m of new accommodation space at 37.1. In the Discussion, we contend that there are significant advantages in this nomenclature system as compared to other systems.

Glacioeustatic Cycles in the Context of Basin Dynamics and Sedimentation

In model/data comparison in Part 4, we will see differences between predicted and observed ages. Here, we use some simple PFMs to demonstrate that net subsidence and sedimentary facies may account for these discrepancies over periods of about one million years.

Figure 12 presents sea-level calculation B for the lower Aptian, with additional graphics representing simple PFMs for basin subsidence and shallow-water carbonate sedimentation (red line) and basin subsidence and deeper-water shale/shallow-water carbonate sedimentation (green line). The diagram thus presents some of the complexity that must be taken into account when dividing stratigraphic sections into various order-cycles and when correlating cycles between stratigraphic sections. Figure 12 should be read in conjunction with Figure 13 that shows the effects of sea-level changes in various parts of a sedimentary basin.

Consider first the sea-level file Figure 12. At this scale, note the complexity introduced by 20-ky (precession), 40-ky (tilt) and 100-ky (eccentricity) cyclicity superimposed upon the higher-amplitude 400-ky cyclicity. Precession and tilt are near ubiquitous as low-amplitude, high-frequency flutter. In contrast, the 100-ky signal sometimes stands out and sometimes blends in with the high-frequency flutter. For the 400-ky peak centered on 121 Ma, 100-ky peaks divide this into two sedimentation events followed rapidly by subaerial exposure and a diagenetic overprint. Elsewhere, such as the 400-ky peaks centered at 117.7 Ma and 122.6 Ma, only precession and tilt high-frequency flutter is discernable. No 100-ky signal is apparent. Thus, the problem of recognizing the 100-ky signal and higher order cycles does not necessarily lie with the sequence stratigrapher. It should be noted that this problem is not affected by increased sedimentation rate. The lack of a consistently discernable 100-ky signal is an attribute of sea-level calculation B.

For the shallow-water PFM (red line in Figure 12), shallow-water carbonate sedimentation is considered as filling new accommodation space during each highstand of the sea. For the deeper-water PFM (green line), the relatively slow deep-water (shale) sedimentation rate increases by a factor of three as shallow-water carbonate sedimentation taken place in water depths of less than 5 m.

Consider the complexity that is presented by the red line with regard to subdivision into various order cycles. Beginning at the 123-Ma highstand, an 800-ky gap separates the previous cycle from what follows. However, in the absence of an input sea-level file, the several 1- to 3-m shallow-water carbonate units that show evidence of subaerial exposure (122.2 to 121.5 Ma), would probably be grouped with the previous Highstand Systems Tract (HST). This in itself has an error of 1.5 my. The unconformity that marks the larger time gap and thus the sequence boundary is at 123 Ma, not 121.5 Ma. Returning to the last subaerial surface of the classical HST (red, 121.5 Ma), there follows the Transgressive Systems Tract (TST) of the next third-order cycle. At 121 Ma, there is a Maximum Flooding Surface (MFS). Everything above the MFS is HST for about 2.5 my. However, without benefit of sea-level input file, the sequence stratigrapher would probably not realize that two 400-ky peaks in sea level did not quite make it back to the subsidence curve (120.5 Ma and 119.7 Ma). Thus, a ‘Fischer plot’ (Fischer, 1964) would have sea level rising whereas the sea-level file indicates it is falling. A major subaerial exposure surface at 118.5 Ma was flooded by a sea-level rise at 117.3 Ma.

In the case of the deeper-water sediments (green line), shallow-water sedimentation occurs at 122.3 Ma and there is subaerial exposure at 121.6 Ma. A HST that culminates in water depth of 2 m at about 119.5 Ma follows the MFS at 121 Ma. A similar sequence of MFS and HST culminates in subaerial exposure at 117.5 Ma.

It is interesting to correlate these two theoretical sections. The shallow-water section consists of one third-order cycle bounded by subaerial exposure surfaces at 121.5 Ma (discussed above) and 118.5 Ma. In contrast, the deeper-water section consists of two third-order cycles (122.3 to 119.5 Ma and 119.5 to 117.5 Ma). Without the benefit of a sea-level file, this would be difficult to resolve. Very likely, the sequence stratigrapher would correlate the subaerial exposure surfaces red 118.5 Ma to green 117.6 Ma, and red 121.5 Ma to green 121.6 Ma. In doing so, he would be miscorrelating the two sections by more than 1 my at the top the sections. With the aid of a sea-level file, the solution of this problem is both simple and intellectually pleasing.

The number and temporal distribution of sedimentation events in shallow water is governed to a great extent by subsidence rate. In Figure 12, if the subsidence rate were less, that is the slope of the red line were a little less, shallow-water units of age 122.2 to 121.1 Ma might not be deposited. In this case, correlation of subaerial exposure surfaces (now red 123 Ma to green 121.6 Ma) would be miscorrelating the two sections by 1.4 my at the base of the sections. Elsewhere, depending on subsidence rate, 400-ky depositional units above the 121.0-Ma highstand might number anywhere from eight to two. Similarly, given a slightly smaller subsidence rate, the deep-water section might transgress the basin margin and become a shallow-water section at the 119.5-Ma intermediate lowstand. This is an example of the utility of a credible high-frequency sea-level file as compared to an ad hoc empirical compilation of past experience.

In Part 4 and Discussion below, most things will be stated as though shallow-water sedimentation fills accommodation space during each peak of the sea-level file. The example provided in Figures 12 and 13 should be regarded as a universal caveat.


We now compare our predicted sequence stratigraphy (columns 3 to 7 in each panel of the Enclosure) to observational data summarized in Sharland et al. (2001) (column 8) and de Graciansky et al. (1998) (column 9). Sharland et al. emphasized the recognition of Maximum Flooding Surfaces (MFS). Their observations are directly comparable to our color-coded picks of new accommodation space in columns 6 and 7. Papers in de Graciansky et al. (1998) focused mostly on identification of stage boundaries and Sequence Boundaries (SB), both of which we regarded as sea-level lowstands (blue in column 3 represent third-order lowstands and blue dots in column 4 represent fourth-order lowstands).

For most of the MFS under consideration here, Sharland et al. (2001, Table 2.3) listed the ‘main driver’ as eustasy. Seven MFS were listed in Table 2.3 as “subsidence”, but the more expanded description in Chapter 4 listed all seven as, “subsidence enhanced by eustasy”. We agree with this definition of the driving mechanism. In the section on Model/Data Comparison below, we address the question of whether the driving mechanism is more explicitly, “ … enhanced by glacioeustasy. ”

A special note should be made of the ages of the various MFS in column 8 as assigned by Sharland et al. (2001). They are not absolute but best guesses from within the defining biozone, and the reader is referred to pages 26 and 51 of Sharland et al. (2001) for a discussion on their uncertainty. In column 8, dates are given only to the nearest million years and, as a general rule, the accuracy (+) assigned must be at least as large as that given to the nearby stage boundaries (column 9). However, the situation is complex and many MFS bring offshore fauna into the local section only briefly. The paleontologist gets a glimpse of a fauna and can pick a biostratigraphic zone but exactly where the fauna is within the zone may be unclear. This is a question of precision within the context of accuracy of the geologic time scale (column 9). We do not fault the use of MFSs (e.g. Sharland et al. (2001) but we suggest that both MFS and SB should be considered. Recognizing MFS is a good way of subdividing a sequence, but their ages are often uncertain. In contrast, the ages of SB are less subject to uncertainty.

When comparing our model results to stage boundaries (column 9), note the caveat regarding physical stratigraphic and biostratigraphic definitions of the boundaries. In the sense of physical stratigraphy, we seek to identify those glacioeustatic lowstands of sea level that might produce unconformity-bounded classic type sections for which the stages were named. Although biostratigraphers begin with these, they often attach boundaries to first and last appearances of taxa in deeper-water sections, far removed from the type sections. Most of the papers in de Graciansky et al. (1998) displayed both biostratigraphic and physical stratigraphic definitions. We color-coded stage names (column 9 in each panel of the Enclosure) with regard to the fit of our predictions to the observational data: green, good; yellow, fair; red, poor.

When considering each time interval, we first noted the major characteristics of the model prediction. Next, we inspected for systematic differences between model and observed time scales. Points correlated in this comparison are shown in bold face and underlined in columns 5 and 9. We then compared model and data regarding MFS and SB.

Model/Data Comparison

We recognize that parameters other than glacioeustasy can cause stratigraphic discontinuities and cycles over a wide range of observational time scales (see Figure 1), but we do not consider them in the context of this paper. The question here is model/data comparison, not whether our model pre-empts all other possible models. We begin at the top of the section because the precision of both model and data is expected to decrease with increasing age. Of particular interest are the Sequence Boundaries. We use some of these as calibration tie points (Ma5, Sa1, Tu4, and so on; see Enclosure column 9) for comparing the orbital-forcing model with the observational data. On the accompanying Enclosure, a double line links the observed Sequence Boundary with its corresponding pedicted lowstand.

Maastrichtian, Campanian and Santonian

In columns 4 to 7 of Enclosure Panel 1, note the following:

  • third-order lowstands that are prominent at 29.0 and 36.0 fourth-order hemi-cycle ID numbers;

  • potentially recognizable third-order highstand signatures for cycles 30, 31, 33, 35 and individual spike highstands at 32.7 and 34.5;

  • Campanian and Santonian (third-order cycles 31-36) have a generally low-amplitude sea-level signature (standard deviation of 3 m) when compared with Maastrichtian above and Coniacian and older below.

Time-scale comparisons: The next step in comparing prediction to observation was to investigate the data for systematic differences in time scales. We correlated Sequence Boundary tie point Ma 5 with our 29.0 and thereby found that the observational time scale was about 200 ky older. Likewise, Sa 1 was correlated with our 36.0 and showed that the observational data is about 500-ky older than our orbital forcing time scale. Thus, when comparing orbital forcing predictions to observations of Maximum Flooding Surfaces and Sequence Boundaries we can expect the data in columns 8 and 9 of the Enclosure to be correlated with events slightly younger in columns 4 to 7.

Maximum Flooding Surfaces (column 8): MFS K180, K170, and K160 are interpreted as correlating with initial flooding events 29.1, 33.1, and 36.1. Sequence Boundaries Ma1 to Ma4 appear to be equivalent to all third- and fourth-order cycles predicted for this interval, and five fourth-order sequence boundaries were predicted and four are observed. We correlated Cam9 with our 31.0; Cam8 with 32.0; Cam6 with 33.0; Cam4 with 34.0; Cam1 with 35.0. Sa1 equates to our 36.0 and we considered that the Santonian consists of a single third-order cycle. Sequence boundaries not mentioned above are thought to represent intermittent detection of the 400-ky signal. Observations picked up 70 percent of predictions.

Coniacian, Turonian and Cenomanian

In columns 4 to 7 of the Enclosure, note the following:

  • third-order lowstands that are prominent at 37.0, 38.0, and 42.0;

  • several fourth-order lowstands that represent persistent returns to maximum ice volume from 40.2 to 41.0;

  • highstands throughout the interval are low compared to older and younger intervals;

  • potentially recognizable third-order highstand signatures in cycles 37, 39, and 40 and an individual fourth-order spike highstand at 41.3;

  • Coniacian through Cenomanian third-order cycles 37 to 42 have generally higher-amplitude sea-level signatures (standard deviation 6.5 m) than the Campanian and Santonian (standard deviation 3 m).

Time-scale comparisons: Tu4 was correlated with our 38.0 and indicated that the observational time scale was about 800-ky older at this tie point. We were particularly interested in the Turonian/Cenomanian boundary because Gradstein et al. (1994, 1995) claimed higher precision here than for boundaries above and below. We correlated Ce5 with our 40.0 and found the observational time scale to be about 300-ky older at this tie point. With regard to the coincidence of our 41.3 with K130 and Ce4, our prediction was that third-order sequence 41 would be characterized by unusually high sea level at 41.3, greater than anything for 1 my on either side. We therefore correlated Sequnce Boundary Ce4 with our 41.2 lowstand. We considered our prediction to be consistent with the observations and cite this as an exact tie point between time scales. Therefore, when comparing orbital forcing predictions with observations of MFS and SB the data in columns 8 and 9 of the Enclosure are correlated with events slightly younger in columns 4 to 7 above our 41.2, whereas direct correlations occur below 41.2 until further notice.

Maximum Flooding Surfaces (column 8): We correlated MFS K150 with our post-Turonian flooding events of early to middle sequence 38. We interpreted K140 as being consistent with our progressively deepening third-order signature and correlated it with our 40.5 and K130 with 41.3. We tentatively correlated K120 with 42.1 or 42.3, but we note the disagreement within published literature regarding the age of Ce1. We therefore have no opinion at this stage regarding the Cenomanian/Albian boundary.

Sequence Boundaries (column 9): The recognition of a single SB within the Coniacian (Co1) was consistent with our third-order sequences 38 and 37. We correlated Tu4 with our 38.0 and Tu3 and Tu2 appear to be related to fourth-order cycles that correlated well with our highstand events 39.7 and 39.5. We correlated Tu1 with our 39.0, Ce5 with 40.0, and Ce4 with 41.2. Ce3 appeared to have a reasonable fit with our 41.0 and we tentatively correlated Ce1 with our 42.0. We note that authors in de Graciansky et al. (1998) placed Ce1 at 99 Ma. However, we note the age of 98.5 Ma assigned to our 42.0 as being the average of the two competing estimates. We tentatively promote our 42.0 as being consistent with the prevailing definitions of the Cenomanian/Albian boundary and suggest that there may be some questions still unanswered regarding the sequence stratigraphic definition of this boundary (see below).

Albian, Aptian and Barremian

In columns 4 to 7 of the Enclosure, note the following:

  • fourth-order lowstands from 43.0 to 46.12 represent persistent returns to maximum ice volume;

  • lower highstands throughout the interval as compared to older and younger intervals;

  • fourth-order lowstands from 47.4 to 48.8 that are substantially lower than lowstands of surrounding intervals, but do not reach levels representing maximum ice volume;

  • highstands throughout this interval are low as compared to older and younger intervals;

  • single prominent lowstands at 48.0, 49.0, 50.8, 50.0, 52.4, 52.0, 54.4, and 54.0;

  • potentially recognizable third-order highstand signatures for cycles 43 and 46 to 54, and an individual fourth-order spike highstand at 44.3.

Time-scale comparisons: As noted above, we correlated Ce4 with our 41.2 and considered the two time scales to be in agreement at this tie point. For reasons explained below, we believe that the position of the Cenomanian/Albian boundary is not yet settled and we will not persue the matter further. Next, we considered Al2 and Ap6 to be synonymous and to correlate with our 48.0. An age of 112.6 + 0.2 Ma encompasses these three age estimates. We correlated Ba0 with our 54.0 and cite this as a tie point for the time scales. Thus, we expect all correlations for this time interval to be horizontal lines.

Maximum Flooding Surfaces (column 8): Comparison of our Albian results with those of Sharland et al. (2001, fig. 4.59) and Immenhauser et al. (1999) indicated several problems. Firstly, the Sharland et al. definition of K100 (Table 4.95) referenced Marker Limestone beds II, whereas Figure 4.59 has K100 labeled as the SB8/9 MFS of Immenhauser et al. (1999). Sharland et al. (2002, Errata, GeoArabia, v. 7, no. 1, p. 127) confirmed K100 is the SB6/7 MFS of Immenhauser et al. (1999). Thus, we tentatively identified K110 as correlative with Immenhauser et al. (1999) SB8/9 MFS and our 44.3 and K100 with 45.5, but we question the age assignments based on the discrepancy noted above. As depicted in Figure 4.59 of Sharland et al. (2001), we correlated K90 with our 47.9 and 47.7 that had been depicted by Immenhauser et al. (1999) as two higher-order flooding events occurring close together. There remains the question of three additional MFSs identified by Immenhauser et al. (1999) in this section. They are shown in Figure 4.59 of Sharland et al. (2001) but not labeled. We comment further on this in the Regional and Local Model/Data Comparisons section.

With regard to Aptian/Barremian MFS, we correlated K80 with our 50.5; K70 with 52.11 to 52.7; K60 with 53.7-53.5; and K50 with 54.5. Since the ages of these MFS are only given to the nearest million years, we consider K80 through K50 to be in good agreement with our sea-level predictions.

Sequence Boundaries (column 9): We were especially encouraged by the fact that sequence boundaries exceeded the number of third-order cycles for this interval. This implies the stratigraphic data are picking up numerous fourth-order cycles. Especially good agreement occurs in lower Aptian and lower Barremian. We correlated Ap3 through Ba6 with our 51.0 to 52.0 and Ba5 through Ba0 with 53.0 to 54.0.

With regard to stage boundaries, we see a particular problem in attempting to reconcile our sea-level file with the common positioning of the Albian/Cenomanian boundary. As these type sections were originally defined as unconformity-bounded, we find it hard to believe that the flooding following our event 42.0 would have been identified but not the flooding that followed event 43.0. In favor of the data from our sea-level file, the last event in the Albian of Oman is indeed flooding and the appearance of planktonic foraminifera. But why is the flooding following our 43.0 not as well known as that following 42.0? We shall return to this under Summary of Time-Scale Comparisons below.

Hauterivian, Valangenian and Berriasian

In columns 4 to 7 of the Enclosure, note the following:

  • third-order lowstands that are prominent at 55.0, 56.0, 57.0, and 61.0;

  • several fourth-order lowstands occur from 58.0 to 60.10;

  • most highstands in this interval are low, as compared to older and younger intervals;

  • potentially recognizable third-order highstand signatures for cycles 55 through 60/61.

Time-scale comparisons: As noted above, we correlated Ba0 with 54.0 and cited it as an exact tie between time scales. Next, on the basis of the ‘second-order’ transgressive/regressive relationships of Jacquin, Rusciadelli et al. (1998, Figure 1), we correlated Ha2 with our 56.0 and find the observational data to be only about 200-ky older at this tie point. This also is an excellent result given that the age uncertainty in this time interval is + 1.9 my. For reasons described below, we correlated Be2 with our 61.0 and found the observational time scale to be only about 410-ky younger at this tie point. Thus, we expect all correlations for this time interval to be horizontal lines.

Maximum Flooding Surfaces (column 8): We correlated MFS K40 with our 55.1, K30 with 58.3 and/or 58.1, K20 with 59.3 and/or 59.1, and K10 with 61.5.

Sequence Boundaries (column 9): Observational data again identified more sequence boundaries than our third-order cycles, and the numbers increased with age. The observational data picked up three sequence boundaries on average, whereas we predicted that each third-order cycle should consist of three to five recognizable flooding surfaces separated by fourth-order sequence boundaries. We correlated Ha4 with our 55.0, Va4 with 57.0, and Be7 with 59.0. We deal with the Be/Ti (Jurassic/Cretaceous) boundary below and for comments on other stage boundaries see the discussion of the lack of systematic differences between time scales above.

Tithonian and Kimmeridgian

Evidence is mounting that a very large meteorite impact event at the Jurassic/Cretaceous boundary produced the Morokweng crater in South Africa (Koeberl et al., 1997; Montanari et al., 1998). When we correlated our cycle 62 with the observational interval Ti5 to Be1, the duration was 2.02 my predicted, versus 800 ky observed. This was the largest discrepancy in our study by a factor of two. We ascribed truncation of cycle 62 to meteorite impact during our 62.4, and correlated Be1 with that event.

In columns 4 to 7 of the Enclosure note the following:

  • third-order lowstands at 62.0, 65.0, and 66.0;

  • several persistently low fourth-order lowstands from 63.0 to 64.8;

  • highstands are low in this interval as compared to those of older and younger intervals;

  • potentially recognizable third-order highstand signatures occur in virtually all cycles from 62 to 66.

Time-scale comparisons: As noted above, we correlated Be1 with our 62.4 and found that the observational data was about 400-ky younger than our orbital forcing time scale. Similarly, Ox 8 was correlated with 66.0 and was found to be about 500-ky younger. These differences are small relative to the + 3-Ma age uncertainty for this time interval. Thus, in comparing orbital forcing predictions to observations of maximum flooding surfaces and sequence boundaries, we will expect all correlations for this time interval to be time equivalent.

Maximum Flooding Surfaces (column 8): MFS J110 through J60 present a special problem as we predicted seven flooding events from 64.7 to 65.1, whereas Sharland et al. (2001) recognized only four. The four fourth-order maximum flooding surfaces observed were assigned ages that concentrated them within our third-order cycle 65. In contrast, we would expect the signature of our cycle 64 to be well represented in the stratigraphy of the Arabian Gulf region. We considered our third-order signatures (column 4) from 62.5 to 66.1 to be correlative with Saudi Arabian formations Sulaiy, Hith, Arab A-D, and Jubaila (Al-Husseini, 1997). We will discuss this further under Regional and Local Model/Data Comparisons below.

Sequence Boundaries (column 9): We correlated Ti5 with our 62.0, Ti2 with 63.0, and Ti1 with 64.0. Although Jacquin, Dardeau et al. (1998) did not provide magnitude information for each third-order cycle in this interval, they did provide ‘second-order ‘ regressive/transgressive information on Tethyan stratigraphic cycles. It is of interest that they characterized everything below Ti 5 as being regressive and everything above as transgressive. This is in excellent agreement with our major flooding at 62.1 and our generally declining sea level for highstands 63.1 to 63.9. Sequence boundaries Ti1 through Ki1 appeared to be intermittently picking up fourth-order cyclicity correlative with our 64.6 through 66.2 (see also Regional and Local Model/Data Comparisons below).

Oxfordian, Callovian and Bathonian

In columns 4 to 7 of the Enclosure note the following:

  • strings of fourth-order lowstands represent near-maximum ice volumes at 66.0 to 69.12, 70.0 to 71.10, and 72.0 to 73.10;

  • as elsewhere, highstands throughout these intervals are relatively low as compared to those of older or younger intervals;

  • prominent single, third-order lowstands at 69.0 and 71.0;

  • potentially recognizable third-order highstand signatures for cycles 68, 69/70, and 71/72.

Time-scale comparisons: As noted above, we correlated Ox8 with our 66.0 and found the observational data to be 500 ky younger than our orbital forcing time scale. Similarly, Ox1 and Ox0 were correlated with 68.0 and 69.12 and found to be about 200-ky younger at this tie point. Next, we correlated Ca 0 with 70.0 or 71.12 and found the two time scales to match with an uncertainty of + 400 ky; and Bt 1 and Bj 5 were correlated with 72.0, 73.12, and/or 73.10 at an accuracy of + 500 ky. Thus, in comparing orbital forcing predictions with observations of maximum flooding surfaces and sequence boundaries, we will expect all correlations at this time interval to be time equivalent.

Maximum Flooding Surfaces (column 8): We correlated MFS J50 with our 67.1 or younger. On the basis of Sharland et al. (2001, fig. 4.41), we correlated J40 with our 69.7 or with our 70.3. In either case we have an age discrepancy of about 2.0 my. J30 was correlated with our 72.3.

Sequence Boundaries (column 9): Ox5 through Ox0, made excellent correlations lowstand for lowstand with our 67.0 through 69.12 (excluding minor lowstand 68.4). Sequence boundaries Ca5 through Bt1 appeared to be related intermittently with fourth-order cycles, but we did not see definite comparisons except in the case of stage boundaries (see comparison of time scales above).

Bajocian, Aalenian and Toarcian

In columns 4 to 7 of the Enclosure we see new patterns:

  • each third-order cycle (with the exception of 73) begins with a prominent lowstand.

  • potentially recognizable third-order highstand signatures occur in cycles 73 to 81.

Time-scale comparisons: The Bathonian/Bajocian boundary was defined with an accuracy of + 500 ky, as noted above. By correlating Bj1 with our 75.0 the observational time scale was shown to be in agreement about + 200 ky and To 7 correlated with 77.0 and indicated that the observational data was about 500-ky younger at this tie point. On correlating the Pliensbachian/Toarcian boundary with 81.0, the two time scales were found to be within about + 150 ky. These time-scale differences are within the age uncertainty for the time interval; thus, when comparing orbital-forcing predictions to observations of maximum flooding surfaces and sequence boundaries, we can expect all correlations for the Bajocian to Toarcian to be time equivalent.

Maximum Flooding Surfaces (column 8): We correlated MFS J20 with our 75.5 and noted that in Sharland et al. (2001, Figure 4.38) there are approximately the correct number of similar organic-rich shales to match our predicted flooding surfaces 75.5 to 78.1. Similarly, we correlated J10 with 79.5.

Sequence Boundaries (column 9): We correlated Bj4, Bj3, Bj1, Aa2, To7, To6 to our third-order lowstands 73.0, 74.0, 75.0, 76.0, 77.0, and 78.0. This left Bj2, Aa1, and To5 to be correlated with miscellaneous fourth-order lowstands nearby.

Summary of Time-Scale Comparisons

Statistical comparisons between orbital forcing and observational time scales are important for two major reasons. First, there is the question of numerology versus science. We predicted fourth-order cycles with periods of 400 ky and compared them with observational data that had an uncertainty of 500 ky in the Late Cretaceous but as much as 3 to 4 Ma in the Jurassic. However, we ran the risk of having a prediction for every event, whatever it’s actual cause. We circumvented this problem by first analyzing the duration of third-order cycles and then analyzing precision and lead/lag relationships regarding stage boundaries.

The second reason for giving serious attention to time scale comparisons was our wish to see fruitful interaction between improvements to geological time scales and future astronomical calculations. There is opportunity for synergism here. Better calibration of astronomical calculations at coarse time intervals will allow them to guide geological interpretation at smaller time intervals. Given systematic discrepancies, a work plan for geologists is self-evident. Does some lack of precision lie simply with the sample intervals of biostratigraphic observations? Passing that hurdle, it is better to have high sampling-density records across boundaries at ten carefully chosen, well-studied locations, rather than only one or two. We see this as the way forward.

Duration of Third-order Cycles

As a test of science versus numerology, we compared the predicted duration of third-order cycles (blue in column 3 of Enclosure) with sequence boundaries (e.g. Ma 5 to Be 2, column 9). We excluded from this comparison parts of the Albian and Kimmeridgian as being problem intervals discussed below, and we also excluded those predicted third-order cycle boundaries for which there are numerous fourth-order cycle boundaries. Finally, we deleted the Jurassic/Cretaceous boundary from this analysis because of a probable meteorite impact effect. We were left with 32 prediction/observation pairs for an analysis of the duration of third-order cycles.

Most of our predicted values were 2.0 my or 2.8 my, being one fourth-order cycle to either side of the 2.4-my (g4-g3) modulation (Figure 3). Two predictions went one fourth-order cycle longer than 2.8-my and two predictions were one fourth-order cycle shorter than 2.0 my. The observational data agreed with two of these predictions. Overall, there was general agreement that some third-order cycles are longer and some shorter. The correlation coefficient between predicted and observed third-order cycle duration was 0.64. Given 32 observations, a relationship between predicted and observed existed at a greater than 99 percent confidence level (Dixon and Massey, 1957, Table A-30a). A correlation coefficient of the order of 0.3 to 0.8 would be expected for an infinite population (Dixon and Massey, 1957, Table A-27).

With regard to age differences at the third-order cycle boundary level, 39 pairs of predicted/observed qualified for this analysis. Over the entire study, the average observed third-order cycle boundary was 16-ky later than the model prediction.

Stage Boundary Comparisons

Next, we inspected for lead/lag relationships at stage boundaries. The key assumption here was that the definition of stage boundaries embodies the difference between fourth-order cycle boundaries and the more important third-order cycle boundaries. Not all of our predicted third-order boundaries have been named stage boundaries. However, for this analysis we specified that every named stage boundary would be correlated with one of our predicted third-order cycle boundaries (blue in column 3 of Enclosure).

Working our way down section, the observed ages of Cam/Ma through Co/Sa are consistently older than predicted, averaging about 1 Ma. Tu/Co and Ce/Tu are much more nearly in agreement. Al/Ce (as discussed below) through Ha/Ba observational data are consistently younger than predicted, averaging about -0.96 Ma. Va/Ha and Be/Va are older than predicted by approximately 0.75 Ma. Ti/Be and Ki/Ti are one 400-ky cycle younger than/older than predicted. Ox/Ki is 1 Ma younger than predicted. Ca/Ox through Aa/Bj are consistently older than predicted by about 300 ky. The To/Aa and Pl/To boundary ages are younger than predicted by about 500 ky.

The most striking systematic relationship within this analysis is the shift in general tendency from observed older than model (Sa/Co and above) to observed younger than model (Al/Ce and below).

As noted above, the Cenomanian/Albian boundary poses a special problem. Why have previous workers chosen to place this Sequence Boundary at our inconspicuous 42.0 (98.5 Ma) when our 43.0 (100.5 Ma) SB and 43.3 (100 Ma) MFS appear so appealing? Is our transition 43.0 to 43.3 indeed as dramatic as it appears in sea-level calculation B?

To investigate this problem, we undertook sensitivity tests that allowed greater amplitude to the ice volume. We concluded that our dramatic high sea levels at 43.1 to 43.9 are robust. As defined by unconformable contacts in shallow-water sequences, we stand by 43.0 as the proper position of the Cenomanian/Albian Sequence Boundary.

Interestingly, van Buchem et al. (2002; especially Figures 6, 7 and 14) provided Cenomanian stratigraphy of the Natih Formation of Oman that is relevant to this problem. The Khatah-1 well contains organic-rich intrashelf basinal deposits at about 1,600 m with a gamma-log signature that in general is stronger than elsewhere in the section. This is their third-order MFS I-2. They correlated it with K120 of Sharland et al. (2001) but assign it a “latest Albian (?)’ age, whereas Sharland et al. assigned it to the “earliest Cenomanian.” The discussion here is about the precise age of a biostratigraphic event recorded in deep-water sediments where sedimentation was continuous. We contend that MFS I-2 equals K120 equals our 43.3, and the related SB equals our 43.0 to comprise the Cenomanian/Albian defining eustatic event within shallow-water sequences. Thus defined, the ‘event’ spans about 0.6 my. The precise timing of the biological first or last appearance event elsewhere in the world is irrelevant here. With the Cenomanian thus defined, we predicted at least 16 fourth-order cycles whereas van Buchem et al. (2002, Figure 14) recognized 8. Assuming a subsidence rate of approximately 18 m/my, only one 400-ky cycle is eliminated as overwhelmed by subsidence. When we try to pick ‘flooding surfaces’, we get seven whereas eight were observed. However, the difference between a fourth-order maximum flooding surface and any other fourth-order highstand is a matter of only a few meters. Persistent shallow-water facies (e.g. Salakh-1 section, 100-230 m) lend credence to our persistent excursions to maximum ice volume, 40.2 to 41.0. We consider this a good first pass at model/data comparison. We have correlated the van Buchem et al. (2002) third-order Sequence III/IV unconformity with our 40.0, third-order MFS II-3 with 41.3, and third-order MFS I-2 with 43.3.

Regional and Local Model/Data Comparisons

Nahr Umr Formation: Albian of Oman

This formation is being studied both in the field and in collaborative model/data comparisons. Subaerial exposure surfaces and third- and fourth-order maximum flooding surfaces (both MFS and mfs) abound. The relation of maximum flooding surfaces to subaerial exposure surfaces is especially interesting—see for example Immenhauser et al. (1999, Figure 7) also reproduced as Sharland et al. (2001; Figure 4.59). In the older part of the section, maximum flooding surfaces tend to occur in the middle of sequences but above SB5 (Figure 4.59 of Sharland et al., they are more prevalent immediately above subaerial exposure surfaces. This is the relationship to be expected if subaerial exposure surfaces occur at lowstands early in the section and at highstands later. The current working hypothesis is that SB5 represents the profound sea-level shift between our 46.5 and 46.12. Confirmation of this event lends credibility to sea-level calculation B.

Shu’aiba Formation: middle Cretaceous of Abu Dhabi

Fischer et al. (1995) provided an excellent example of seismic data across the Bu Hasa and Asab Bab fields. Figures 4.57a,b of Sharland et al. (2001) also dealt with these rocks. Carbonate build-up appears to have begun at about MFS K70 (120 Ma) and to have been completed by K90 (111 Ma). There are at least 23 fourth-order cycles organized into five third-order cycles in this time interval. Fischer et al. (1995) appeared to have divided the section into six units, each consisting of a shallow-water and a deeper-water facies; these are most likely third-order cycles. The transition from deep to shoaling Bab is likely to have occurred at lowstand, as indicated in Figure 12. The implied complexity of fourth-and higher-order sedimentation and related early vadose and phreatic diagenesis remains largely unresolved.

Hith and Arab Formations: Upper Jurassic of Abu Dhabi and Saudi Arabia

The Arab Formation deserves special consideration because of its economic importance. Additionally, this is the part of the stratigraphic record in which Sharland et al. (2001) identified the most closely space MFSs. We strongly agree with them that many third- and fourth-order MFS should be recognizable for this time interval in these facies. However, we disagree on some age assignments and the relationships of these named MFS to the stratigraphy of the Arab and Hith formations.

Part of the problem resides with the definition of MFS J110 and J70 in Yemen and in the correlation from Yemen to the Arabian Gulf. Sharland et al. (2001, Figure 4.45) depicts Yemen marls and limestones that contain numerous ammonite-bearing limestones, presumably related to clear-water sedimentation associated with maximum flooding in a shelf section otherwise influenced by land-derived fine-grained clastics. Nine ammonite-bearing limestones are depicted but only three are given MFS number designations (J70-J90). Indeed, the designation J90 is given to a group of three such limestones, suggesting that identification of fourth-order flooding surfaces is possible within these sections. We accept these Yemen strata as an excellent recording device for MFS but we question their correlation to the Arabian Gulf.

With regard to J110, Sharland et al. (2001, Figure 4.43 and Table 4.76) correlated J110 with the uppermost Hith Formation. We agree that an MFS occurs near the top of the Hith but this correlation gives us an age discrepancy of 1.77 Ma. Alternatively, we also predict an MFS near the base of the Hith and this correlation to J110 gives us an age discrepancy of only 450 ky. Thus, we accept the age of J110 as defined by ammonite biostratigraphy in Yemen and correlated it to the base of the Hith Formation. Likewise, we accept the Yemen definition and age of J70, but we will correlate it differently in the Arabian Gulf.

Working down section from the top of the Hith through the Arab Formation to J60 of Sharland et al. (2001) we noted the excellent graphics and descriptions by Azer and Peebles (1998) regarding the Abu Dhabi offshore section. We took the signature of our cycle 63 to be conducive to the accumulation of a thick Hith evaporite. We correlated our 62.1 with the limestone at the top of the Hith Formation and 63.1 with J110 and its base, and, consistent with Sharland et al. (2001, Figure 4.50), we identified five fourth-order cycles within the evaporite portion of the Formation. Next, we took Figure 4.48 of Sharland et al. (2001) and attempted to identify an appropriate number of fourth and higher order cycles proportional to the thicknesses indicated by the variable duration of Arab A, B, and C. We tentatively accepted J100 as the MFS that separates Arab A from Arab B and correlated it with our 64.5, and similarly with J90 (our 64.1) as the MFS that separates Arab B from Arab C, and J80 (our 65.5) as the MFS that separates Arab C from Arab D. This brings us again to J70. We accept the definition and age from Yemen and correlated it with our 65.1.

J60 is rated by Sharland et al. (2001) as a subregional MFS rather than being one of the many plate-wide MFSs. Likewise, they identified correlation confidence for this MFS as only “ moderate” Tables 4.65 and 4.66 give the reference section as the “upper Hanifa shales”. However, Figure 4.43 plots J60 as more like lower Hanifa and Figure 4.44 plots it as the initial MFS of the lowermost Jubailah Formation. We correlated J60 with any of three predicted MFS within cycle 66 and considered it to be within the Jubailah Formation in Saudi Arabia. As indicated in the Enclosure, we consider our 66.1 to be lowermost Kimmeridgian.

With regard to the contact between the Jubailah Formation and overlying Arab D in Saudi Arabia, we suggest that this contact lies somewhere between our 66.7 and 65.1. Near to the contact, deeper-water Jubailah environment changed as sea level fell and shallow-water facies were established. Whereas shallow water facies below this contact would seldom be subaerially exposed, shallowing-upward cycles above the contact would commonly fill the accommodation space and be exposed in Arab Formation carbonates younger than our 65.1.

Finally, we caution that the definition of boundaries within the Arab Formation may vary regionally. Whereas the Arab C, B, and A of Abu Dhabi (e.g. Sharland et al., 2001, Figure 4.48) are composed of stacks of numerous carbonate/evaporite cycles, the Arab Formation of Saudi Arabia’s producing regions is commonly characterized as carbonate/evaporite couplets (e.g. Figure 4.43). Identifying cycles with total duration proportional to the variation in member thicknesses results in different boundary assignments.

Jurassic of Yemen and Oman

Sharland et al. (2001, Figures 4.41 and 4.45) provided excellent examples of ammonite zonation. The ammonites are concentrated in limestone intervals within a shale section but there are many more breaks than there are Sharland-numbered MFS. Likewise, Figure 4.38 provides an excellent example of numerous candidates for additional MFS that have not been named. As with the Albian of Oman (above), many fourth-order maximum flooding surfaces appear to have gone over looked or ignored.


Time Scales

We note that the model/data comparisons made here are no proof of accuracy of the observed or the astronomical time scale. Gradstein et al. (1995) made an excellent estimate of possible errors of boundary ages for the observed time scale and these remain in force. Even though slight changes in the values of the astronomical constant g can produce alternative astronomical time scales that differ from that presented here by millions of years, such discrepancies are typically within the error ascribed to the observed time scale by Gradstein et al. (1995).

Orbital forcing of glacioeustasy is both expected and periodic, and is therefore predictable. In contrast, large meteorite impact events are expected, but not periodic; therefore their timing is not predictable. By examination of Mars and Earth’s moon, it is estimated that 75 to 90 percent of large Mesozoic meteorite impacts on Earth have yet to be found (Grieve, 1997). Technology presented here picks up the South African Morokwong meteorite impact (Cretaceous/Jurassic boundary) (Montanari et al., 1998; Koeberl et al., 1997) as a shorter-than-predicted third-order cycle. Similar shortening of third-order cycles was seen at other times, but was less dramatic. This may provide stratigraphers with time horizons at which they may expect to find evidence of meteorite impact. This is especially important because so much of the Earth is covered by water. Whereas 33 large and well-dated meteorite impact events have been recognized on the continents (Montanari et al., 1998), only one has been found in the ocean basins (Gersonde et al., 1997). Sedimentologic evidence of oceanic meteorite impact is far more likely to be found than the actual crater.

Under Stage Boundary Comparisons (above), we pointed out systematic differences between the astronomical time scale presented here and the observed time scale of Gradstein et al. (1995). We would like to think that astronomers will take up the challenge, but we recognize that few of them are interested in this problem and younger time scales have priority. Thus, results presented here are both a working hypothesis and a target for future investigation.

Eccentricity Modulation at 400 ky and 2.4 my

A graphics problem exists when comparing modulation at these two periods. Modulation looks more important in Figure 3 than in Figure 2. As noted earlier, this is an illusion arising from the longer period of the 2.4-my modulation of the ‘100-ky’ signal. In fact, amplitudes at anti-node are additive and amplitudes at node are subtractive. Thus calculated, (g2-g5) modulation is 29 percent more pronounced than the (g4-g3) modulation. This dictates that large sea-level events can occur within the 400-ky eccentricity cycle when eccentricity modulation interacts with tilt modulation.

Basic to our sea-level calculations is the premise that the Earth tends to accumulate high-latitude continental ice caps if there is insufficient summer heat to melt the ice. However, a high degree of eccentricity near the summer solstice at perihelion creates warm polar summers as does high tilt under the same conditions. Thus, if modulation anti-nodes of eccentricity and tilt coincide, polar summers are warmer than normal. Conversely, if modulation nodes of eccentricity and tilt coincide then polar summers are cooler than normal.

Zachos et al. (2001) demonstrated the latter case at the Oligocene/Miocene boundary by means of the 400-ky (g2-g5) eccentricity modulation and complex tilt modulation. The result was a conspicuous build-up of ice. The event lasted about 150 ky and was more conspicuous than ice build-up events at nearby 2.4-Ma nodes that did not have the tilt effect.

We have inspected our calculations for this effect. Near-coincidence of 400-ky (g2-g5) eccentricity modulation nodes and more complex tilt-modulation nodes is common, occurring approximately every 2 to 4 my, thereby indicating regular ice build-up. However, the more spectacular near-coincidence of 2.4-my (g4-g3) eccentricity modulation nodes and more complex tilt-modulation nodes, which would provide a longer period of low eccentricity and low tilt, does not occur throughout the 190 my of our calculations.

We have generalized the ‘Zachos effect’ to include both coincidence of eccentricity and tilt modulation at nodes (excessive ice build-up) and at anti-nodes (excessive ice melting). Large events either way can occur within a relatively short time. The effect is to give character to the signature of third-order cycles.

Cycle-order Definition and Nomenclature

The tuning fork of geologic time is the 400-ky (g2-g5) cycle. We have numbered these cycles back from the Holocene (highstands, whole numbers; lowstands, half numbers) and designated them as fourth-order cycles. They are organized into third-order cycles by (g4-g3) modulation. Third-order cycles were also numbered from the Holocene. The break between a cycle and the next youngest was taken as the first major transgression of the younger cycle. In most cases, the length of third-order cycle thus defined was 2.4 + 0.4 my (2.0 or 2.8 my). This time interval is in the realm of meteorite impact events. Here, we have found it convenient to recognize the Morokwong event (Jurassic/Cretaceous boundary) as having caused the shortening of an orbitally defined third-order cycle. Likewise, ice build-up Zachos events might become confused with third-order cycle boundaries, but we encountered no such confusion in this study. Thus, we are confident that the orbital forcing definition of third-order cycles will serve us well throughout the Mesozoic.

Higher-order cycles have presented a problem. At first glance, it would seem logical to continue with the expectations of orbital forcing. A large component of eccentricity is the ‘100-ky’ period; tilt, 40 ky; precession, 20 ky. Why not call these fifth-, sixth- and seventh-order cycles? As discussed with regard to Figure 12, the problem becomes the sporadic recognition of the ‘100-ky’ signal. For most Mesozoic purposes, we suggest combining everything above fourth-order as “higher-order cycles.”

With regard to standardized nomenclature, we suggest capital letters for third-order events, such as Maximum Flooding Surface (MFS) and Sequence Boundary (SB), and lower case for the equivalent fourth-order events. Next, is the question of Transgressive System Tract (TST), Highstand System Tract (HST) and Lowstand System Tract (LST). In the simplest terms, anything younger than a SB and older than the next MFS is the TST, and anything between that MFS and the next younger SB is the HST. The LST is recognized as offset basinward from HST in some classic situations. In practice, HST becomes potentially confusing when applied to deep-water sequences (see Figure 12) where the shallowest water (or even subaerial exposure) occurs at the lowstand. We prefer to return to the ‘inclusive terms of ignorance’, Transgression and Regression, when precise causal mechanism or position in the cycle is in doubt. For example, deep-water sediments in Figure 12 are clearly regressive after MFS, but the use of HST would imply culmination of the sequence at highstand, which is not the case.

Note that our system is in contrast to the objective system applied by many sequence stratigraphers, for example, Vail et al. (1997) and Goldhammer et al. (1991). The objective way to handle cycle order is to divide them into order-of-magnitude categories; third-order, 10 my to 1 my; fourth-order, 1 my to 100 ky; fifth order, 100 ky to 10 ky, and so on. Although we take no fault with early workers in being objective, we see several major problems with continuing this scheme. First, if one ever found an 8-my cycle, it would be unusual, new and exciting but it would not be possible to combine it with (g4-g3) glacioeustasy third-order cycles (2.0 to 2.8 my). Second, 100 ky is a poor time interval with regard to glacioeustasy expectations. It combines (g2-g5) 400-ky cycles with much the longer (g4-g3) cycles, and it splits expected eccentricity cycles (95 to 130 ky) into unnatural groups. Thus, we recommend that the traditional definitions be discontinued.

Utility of Simple Parametric Forward Models

The fact that model/data comparison has yielded promising results lends credence to the Parametric Forward Model (PFM) approach to sequence stratigraphy. Although a single forward model is not a unique solution, it demonstrates that the model is capable of convergence with the data. Different parameters may yield similar results, or future parameterizations may better fit the (improved?) data, but the overall structure described here is likely to remain similar.

To achieve rapid, cost-effective closure on tasks such as reservoir characterization and flow simulation, similar PFMs must be constructed for numerous other characteristics of the geological system. Some of these characteristics are global, some regional, some are local; some are easy to model, and some are complex.

The global characters include tectonoeustasy in addition to glacioeustasy (modeled here). As indicated in Figure 1, changing the volume of the global ocean basin involves many processes. As an exercise in dynamic process modeling, this is probably an intractable problem. However, the role of tectonoeustasy in stratigraphy is surely intertwined with that of glacioeustasy. Data for the last 80 my of sea-floor spreading (e.g. Matthews and Frohlich, 1998, Figure 4) suggests that the sea-floor spreading effect is large (hundreds of meters amplitude), slow (about 1% of the glacioeustatic rates calculated here), and unidirectional for long periods. Therefore, for the last 80 my, sea-floor spreading can be taken into account by a very simple PFM of one column of numbers interpolated from the data. Beyond that, where the amplitude of a stratigraphic effect exceeds the effect ascribable to ice volume and local conditions in numerous stratigraphic studies worldwide, then forward modeling of tectonoeustasy, even in the absence of direct observation, deserves consideration. Ice volume cannot exceed the bounds of geological reasonableness. Thus, large eustatic amplitude implied by widespread stratigraphic observations would be a reason to attempt to model tectonoeustasy on older time scales. The physics is well known and time does not turn physics on its head. The grand total of tectonoeustasy as what is required to exceed the limits of glacioeustasy. It follows that one must define major features by observation and connect them with rates of change prescribed by reasonable physics in order to differentiate tectonoeustasy from glacioeustasy.

Important regional characteristics are tectonic setting, subsidence rates, climate and sediment supply, isostasy, and compaction. Important local characteristics include water depth, salinity, diagenesis, and basin geometry. All can be modeled by use of a simple PFM. However, even when each submodel is kept simple, the human mind cannot predict the outcome of the consolidated megamodel resulting from simultaneous consideration of numerous simple PFMs and the modeling effort must be computer-based. For most carbonates, a one-dimensional model will suffice, for example, SedSim (Matthews and Frohlich, 1998; R.K. Matthews, 2000, unpublished R.K. Matthews and Associates Progress Report: Orbital Forcing of Glacioeustasy (4); Carbonate Diagenesis around Water Tables). However, for clastics, at least a two-dimensional model is required, for example, STRATA-various of Frohlich and Matthews (1991).

The Importance of Glacioeustasy

Relative sea level has fluctuated through geologic time due to a variety of processes (Figure 1) and on a wide range of time scales. The size of the fluctuations would have been different from place to place, depending on the process. In general, the observational data suggest that the global sea level has varied by about 200 m above or below its present level. Yet through it all—if we are correct—the glacioeustasy beat goes on at amplitudes of a few meters to tens of meters and at periodicities described above. In this respect, observational papers are now appearing that recognize rapid sea-level events as likely to be of glacioeustatic origin, even in the former supposed ice-free Cretaceous; for example, Immenhauser et al. (1999) and Gale et al. (2002).

As indicated in Figure 14, the coastal plain and continental shelf occupy about 15 percent of today’s Earth. This is a very flat surface. A 3-m rise in sea level can send the shoreline 10 km landward; 30 m, 100 km. For any time in Earth history, we propose that the Earth has had a similar large proportion of its surface in the zone of the coastal plain and continental shelf. Thus, for example, tectonoeustasy may displace the mean sea level by 200 m, and thereby shore lines by several hundred kilometers. But, within the paleogeographic context of relatively slow tectonoeustasy, a relatively rapid 3-m glacioeustatic sea level rise could still send the shore line 10 km landward, and a 30-m rise lead to a 100-km transgression. Thus, variations in the glacioeustatic component of sea level can have a profound effect on the location and type of sedimentation that occurs.


  • A simple Parametric Forward Model (PFM) has yielded promising results in model/data comparison. Orbital forcing affects Antarctic continental ice volume and therefore the volume of water in the oceans; this affects sea level and therefore stage boundaries, sequence boundaries, and maximum flooding surfaces. The result lends credence to the PFM approach to sequence stratigraphy.

  • The orbital forcing 400-ky cycle is the tuning fork of geologic time. This is the orbital forcing (g2-g5) pure sine wave, the highest amplitude term of eccentricity. It is also a prominent modulation of several ‘100 ky’, high-amplitude eccentricity terms and the natural definition of the fourth-order cycle in sequence stratigraphy.

  • Fourth-order cycles are naturally grouped into third-order cycles by orbital forcing (g4-g3) modulation. Because the node in (g4-g3) is broad relative to the fourth-order cycles, the majority of third-order cycles are 2.4 + 0.4 my; that is 2.0 my or 2.8 my. At the 99 percent confidence level, the relation of predicted to observed cycle length is 0.64 (the correlation coefficient for an infinite population of predicted and observed cycles is between 0.3 and 0.8).

  • The geologic time scale and the orbital forcing (astronomical) time scale have separate sources of error that are approximately of the same order of magnitude. In the long term, the goal should be compatibility of the two time scales. Geologists need to concentrate on the causal mechanisms of sequence boundaries, with special attention to distinctions between normal third and fourth-order orbital forcing (periodic and well understood; as set forth here), Zachos events (probably periodic, but not yet well understood), and meteorite impact events (expected and well understood, but chaotic).

  • Signatures are to be expected within third-order sequences. Given the PFM target presented here, sequence stratigraphers should look into them: time scales may change, but the third-order node count and many third-order signatures will remain intact.

  • A time of transition in sequence stratigraphy is approaching. Just as our understanding of geosynclines changed dramatically with the incorporation of plate tectonics into our thinking, our understanding of sequence stratigraphy will change dramatically with the incorporation of orbital forcing of glacioeustasy. Given improved time scales, third-order signatures, and the ever-increasing computing power, the time is approaching for geologists to abandon classic dataintensive sequence stratigraphy in favor of more rapid, cost-effective parametric forward modeling. Although numerous other PFMs will be required (some global, some local), orbital forcing of glacioeustasy is likely to remain a major consideration. We need not reinvent the wheel each time we recommend a wildcat well or study a new field.


We thank Jacques Laskar for his very helpful visit to Brown University in December 2000. We thank (alphabetical) Sadad Al-Husseini, Jack Allan, Eric Barron, Jeff Chen, Al Fischer, Bill Galloway, Tim Herbert, Adrian Immenhauser and Ross Peebles for helpful discussion/review of early versions of this manuscript. We thank GeoArabia ‘s anonymous reviewers for their useful comments, and GeoArabia ‘s editorial staff (David Grainger, Joerg Mattner and Moujahed Al-Husseini) for their assistance in making the paper more reader-friendly. We are especially grateful for their help with Parts 1 to 3 and the contributions of Figures 1 to 4 and 14. The design and drafting of the final graphics was by GulfPetroLink.


Robley K. Matthews is Professor Emeritus of Geological Sciences at Brown University, Rhode Island, USA and a general partner of R.K. Matthews and Associates. He has over 35years experience in carbonate sedimentation and diagenesis and their application to petroleum exploration and reservoir characterization. His current interests center on the use of computer-based parametric forward models in stratigraphic simulation.


Cliff Frohlich is a Senior Research Scientist and Associate Director of the Institute for Geophysics, University of Texas at Austin, and a partner in R.K. Matthews and Associates. He is best known for his research concerning deep-focus earthquakes. However, he also has an abiding interest in computer modeling of various mechanical and dynamic systems, ranging from bowling balls to stratigraphy of sedimentary basins.



In columns 3 to 7: Green, yellow and red signify ad hoc estimates of quality—green, good; yellow, intermediate; red, poor; blue in column 3 indicates beginning of a new third-order cycle; white in column 3 signifies uncertainty.

In columns 6 to 7: Estimates of accommodation space are due to rise in global sea level only; isotsasy and geodesy need to be taken into account and adjustments will vary widely.

In column 8: Notes (1) to (3) relate to MFS occurrences—(1) plate-wide to (3) local. Notes (4) to (8) relate to MFS criteria—(4) ammonite-bearing unit; (5) open-marine indicators; (6) gamma-ray log maxima; (7) marine unit in sabkha-dominated succession; (8) other/complex (see Sharland et al., 2001).

In column 9: Blue dates indicate boundary age estimates (this paper); black dates and sequence boundary numbers are from other authors. White signifies uncertainty.