Realistic Paleobathymetry of the Cenomanian–Turonian (94 Ma) Boundary Global Ocean

: At present, global paleoclimate simulations are prepared with bathtub-like, ﬂat, featureless and steep walled ocean bathymetry, which is neither realistic nor suitable. In this article, we present the ﬁrst enhanced version of a reconstructed paleobathymetry for Cenomanian–Turonian (94 Ma) time in a 0.1 ◦ × 0.1 ◦ resolution, that is both realistic and suitable for use in paleo-climate studies. This reconstruction is an extrapolation of a parameterized modern ocean bathymetry that combines simple geophysical models (standard plate cooling model for the oceanic lithosphere) based on ocean crustal age, global modern oceanic sediment thicknesses, and generalized shelf-slope-rise structures calibrated from a published global relief model of the modern world (ETOPO1) at active and passive continental margins. The base version of this Cenomanian–Turonian paleobathymetry reconstruction is then updated with known submarine large igneous provinces, plateaus, and seamounts to minimize the difference between the reconstructed paleobathymetry and the real bathymetry that once existed.


Introduction
Understanding past climate is critical for the scientific study of ancient, modern, and future climate change. Computer-generated models, also known as general circulation models (GCMs) are often used to study the plausible dynamic effects of diverse boundary conditions on climate. Paleoclimate data provide a useful framework from which these models can be used to predict past, present, and future climate change scenarios. The use of numerical models to investigate Earth's climate history has increased significantly over the past decade [1][2][3]. These models have grown in sophistication, with increased space-time resolution, coupled ocean-atmosphere GCMs that incorporate vegetation, soil (weathering), ice, and chemistry modules.
Accordingly, more robust and detailed boundary conditions are needed. One of the most difficult and often overwhelming tasks in paleoclimate modeling is the assembly of surface boundary conditions for past time periods [2]. These boundary conditions include but are not limited to continental configurations, topography, ocean bathymetry, ocean chemistry, terrestrial vegetation distribution, atmospheric trace gas concentrations, insolation, etc. Many of these boundary conditions have been the subjects of extensive research, and abundant data are available. A notable exception is paleobathymetry, which is often modeled as a bathtub with steep walls or oversimplified bathymetry to such a great extent that the role of bathymetry on climate change is greatly undermined. A simple scholarly search 1.
The occurrence of a major oceanic anoxic event (OAE2) associated with large accumulations of organic rich sediments in the oceans and leading to the largest petroleum resources being utilized by society today [36,37];

2.
Major carbon cycle acceleration and biogeochemical perturbation in the form of a major global marine carbon isotope excursion [38]; 3.
Dramatic global warming from mantle-plume volcanism (main phase Caribbean and Madagascar, late phase High Arctic and Ongtong Java LIPs) [39]; 4.
Extensive mid-ocean ridge formation, high crustal heat flow and rapid ocean basin expansion [15]; 5.
Major change in ocean circulation from eddy-dominant with widespread cosmopolitan biogeography in the Cenomanian to thermohaline with restricted marine biogeographic zones in the Turonian [40]; 6. Highest sea levels of the Mesozoic with development of major, multiple epicontinental seas, some serving as gateways (e.g., the Western Interior Seaway connecting the tropical Atlantic with polar Arctic) [16,41] and 7.
Completion of opening of the Central Atlantic Seaway connecting the North and South Atlantic sectors [42,43]. 8.
Onset of the Alpine orogeny as indicated by the relative plate motions and paleogeographic evidence from Africa, Iberia and Europe (especially, the Carpatho-Balkan and Sicily regions) [44][45][46][47].
All of these events involved the oceans and climate change, and so it is clear that a detailed C-T paleobathymetry model will be an important contribution for future paleoclimatic investigations.

Materials
To reconstruct the C-T bathymetry for the global ocean with generalized continental shelf-slope-rise structures using the OESbathy methodology (version 1.0, [35]) five major datasets were used. A brief description of these individual datasets is provided below.
Ocean crust age: Müller [16] published global reconstructions of ocean crust age in 1 Myr intervals for the past 140 million years. For each reconstructed age in [16], ocean crust age, depth to basement, and bathymetry are given. This reconstructed bathymetry obtained as part of the dataset provided by Müller [16] is referred to hereafter as EB08 (EB: EarthByte). These data are in 0.1 • × 0.1 • resolution (3601 longitude × 1801 latitude points). For this project, (modern) crustal age reconstruction data are used ( Figure 1).
Geosciences 2018, 8, 21 3 of 20 2. Major carbon cycle acceleration and biogeochemical perturbation in the form of a major global marine carbon isotope excursion [38]; 3. Dramatic global warming from mantle-plume volcanism (main phase Caribbean and Madagascar, late phase High Arctic and Ongtong Java LIPs) [39]; 4. Extensive mid-ocean ridge formation, high crustal heat flow and rapid ocean basin expansion [15]; 5. Major change in ocean circulation from eddy-dominant with widespread cosmopolitan biogeography in the Cenomanian to thermohaline with restricted marine biogeographic zones in the Turonian [40]; 6. Highest sea levels of the Mesozoic with development of major, multiple epicontinental seas, some serving as gateways (e.g., the Western Interior Seaway connecting the tropical Atlantic with polar Arctic) [16,41] and 7. Completion of opening of the Central Atlantic Seaway connecting the North and South Atlantic sectors [42,43]. 8. Onset of the Alpine orogeny as indicated by the relative plate motions and paleogeographic evidence from Africa, Iberia and Europe (especially, the Carpatho-Balkan and Sicily regions) [44][45][46][47].
All of these events involved the oceans and climate change, and so it is clear that a detailed C-T paleobathymetry model will be an important contribution for future paleoclimatic investigations.

Materials
To reconstruct the C-T bathymetry for the global ocean with generalized continental shelf-sloperise structures using the OESbathy methodology (version 1.0, [35]) five major datasets were used. A brief description of these individual datasets is provided below.
Ocean crust age: Müller [16] published global reconstructions of ocean crust age in 1 Myr intervals for the past 140 million years. For each reconstructed age in [16], ocean crust age, depth to basement, and bathymetry are given. This reconstructed bathymetry obtained as part of the dataset provided by Müller [16] is referred to hereafter as EB08 (EB: EarthByte). These data are in 0.1° × 0.1° resolution (3601 longitude × 1801 latitude points). For this project, (modern) crustal age reconstruction data are used ( Figure 1).   Three geologic time scales that concern the age of the Cenomanian-Turonian (C-T) boundary are summarized in Table 1 The C-T boundary age is now extremely well constrained at 93.9 Ma [48]. The duration of the interval from the C-T boundary to the top of C34 is estimated to be from 8.875 m.y. (GTS2012) to 8.5 m.y. (GTS2016), i.e., 23% to 24% of the duration of the C34 to M0 interval. This scales the C-T boundary age in the Müller GTS to between 92.375 Ma and 92.0 Ma. This differs slightly from the 93.5 Ma C-T boundary age of Gradstein [49], whose time scale was used to calibrate the digital age grid for the oceanic crust by Müller [31] and later used in EB08 [15]. Our study assumed a "traditional" age of 90 Ma for the C-T boundary for the ocean crustal component (construction of EB08). However, our reconstruction is actually between 2.375 m.y. and 2 m.y. younger than the C-T boundary age according to the Müller GTS (see Table 1), the latter presumably equivalent to the 93.9 Ma age of the C-T boundary assigned by GTS2012 and GTS2016.
Modern ocean sediment thickness: Modern ocean sediment thickness data from Divins [51] and Whittaker [52] are used. These data are derived from seismic profiling of the world's ocean basins and other sources. The reported thicknesses are calculated using seismic velocity profiles that yield minimum thicknesses. Data values represent the distance between the seafloor and "acoustic basement". The data are given in 5s × 5s resolution and have been regridded to 0.1 • × 0.1 • resolution values to match the EB08 grid [35]. ETOPO1: To reconstruct shelf-slope-rise structures, large igneous provinces and seamounts, modern topography and bathymetry from ETOPO1 [53] is used. The "Bedrock" version of ETOPO1 is utilized for this work, which is available in a 1s × 1s resolution (earthmodels.org), and regridded to 0.1 • × 0.1 • resolution [35] in order to match the EB08 grid. This version of ETOPO1 includes relief of Earth's surface depicting the bedrock underneath the ice sheets. However, only the oceanic points in this dataset are used, so that there is no impact on the reconstructed bathymetry.
Large igneous provinces: In this reconstruction, the large igneous provinces (LIPs) dataset [54] is used. This dataset includes a geographic description of each LIP (Figure 2), plate identification number, and times of appearance and disappearance (in Ma). These LIPs are distinguished from mid-ocean ridge and subduction-related rocks based on petrology, geochemistry, geochronology, physical volcanology, and geophysical data. For this study only oceanic LIPs are considered. The modern distribution of the oceanic LIPs ( Figure 2a) were rotated back to 94 Ma (Figure 2b).
Global Seamount Database: In this study, the seamounts from the Global Seamount Database [55] are used, which is based on satellite-derived vertical gravity gradient (VGG) data with globally 24,643 potential seamounts that are higher than or equal to 0.1 km, and located away from continental margins. The modern distributions of these seamounts ( Figure 3a) were rotated back to 94 Ma (Figure 3b).

Methods
For this reconstruction, a standard plate cooling model for the oceanic lithosphere is combined with the age of the oceanic crust, global oceanic sediment thicknesses, plus generalized shelf-sloperise structures calibrated at modern active and passive continental margins. As with the modern ocean basins, there are different types of crust: oceanic crust, submerged continental crust, and a transitional zone between the two (e.g., [35]). In this reconstruction, the regions underlain by oceanic crust to which an age has been assigned are termed "open ocean" regions. Here, a negative sign is assigned to signify depths below mean sea level (MSL).
The reconstruction of open ocean bathymetry starts with ocean crust age ( Figure 1). This information is available only at locations where oceanic crust is preserved or has been reconstructed. The ocean depth-to-basement ( Figure 4) is calculated based on a cooling plate with an average depth of −2639.8 m for the depth of mid-oceanic ridge [35].
The addition of sediment and an isostatic correction from sediment loading of the depth-tobasement oceanic crust (e.g., [56]) is critical to complete the bathymetric reconstruction of the open ocean region. A parameterized multilayer sediment cover, called "OES sediment thickness", based on a third-degree polynomial fit between area corrected global sediment thickness data of modern oceans and marginal seas [51,52] and the age of the underlying oceanic crust [35], was isostatically added on top of the depth-to-basement. The sediment loading approach used here is similar to procedures used by Crough [57] and Sykes [58], was adopted using a multicomponent sediment layer with varying sediment densities calculated from a linear extrapolation of published sediment densities [35]. Once the sediment loading was added, the reconstruction of the open ocean regions was considered complete ( Figure 5), and ready for merging with shelf-slope-rise structures, and adding LIPs, plateaus and seamounts, as will be described next.

Methods
For this reconstruction, a standard plate cooling model for the oceanic lithosphere is combined with the age of the oceanic crust, global oceanic sediment thicknesses, plus generalized shelf-sloperise structures calibrated at modern active and passive continental margins. As with the modern ocean basins, there are different types of crust: oceanic crust, submerged continental crust, and a transitional zone between the two (e.g., [35]). In this reconstruction, the regions underlain by oceanic crust to which an age has been assigned are termed "open ocean" regions. Here, a negative sign is assigned to signify depths below mean sea level (MSL).
The reconstruction of open ocean bathymetry starts with ocean crust age ( Figure 1). This information is available only at locations where oceanic crust is preserved or has been reconstructed. The ocean depth-to-basement ( Figure 4) is calculated based on a cooling plate with an average depth of −2639.8 m for the depth of mid-oceanic ridge [35].
The addition of sediment and an isostatic correction from sediment loading of the depth-tobasement oceanic crust (e.g., [56]) is critical to complete the bathymetric reconstruction of the open ocean region. A parameterized multilayer sediment cover, called "OES sediment thickness", based on a third-degree polynomial fit between area corrected global sediment thickness data of modern oceans and marginal seas [51,52] and the age of the underlying oceanic crust [35], was isostatically added on top of the depth-to-basement. The sediment loading approach used here is similar to procedures used by Crough [57] and Sykes [58], was adopted using a multicomponent sediment layer with varying sediment densities calculated from a linear extrapolation of published sediment densities [35]. Once the sediment loading was added, the reconstruction of the open ocean regions was considered complete ( Figure 5), and ready for merging with shelf-slope-rise structures, and adding LIPs, plateaus and seamounts, as will be described next.

Methods
For this reconstruction, a standard plate cooling model for the oceanic lithosphere is combined with the age of the oceanic crust, global oceanic sediment thicknesses, plus generalized shelf-slope-rise structures calibrated at modern active and passive continental margins. As with the modern ocean basins, there are different types of crust: oceanic crust, submerged continental crust, and a transitional zone between the two (e.g., [35]). In this reconstruction, the regions underlain by oceanic crust to which an age has been assigned are termed "open ocean" regions. Here, a negative sign is assigned to signify depths below mean sea level (MSL).
The reconstruction of open ocean bathymetry starts with ocean crust age ( Figure 1). This information is available only at locations where oceanic crust is preserved or has been reconstructed. The ocean depth-to-basement ( Figure 4) is calculated based on a cooling plate with an average depth of −2639.8 m for the depth of mid-oceanic ridge [35].
The addition of sediment and an isostatic correction from sediment loading of the depth-to-basement oceanic crust (e.g., [56]) is critical to complete the bathymetric reconstruction of the open ocean region. A parameterized multilayer sediment cover, called "OES sediment thickness", based on a third-degree polynomial fit between area corrected global sediment thickness data of modern oceans and marginal seas [51,52] and the age of the underlying oceanic crust [35], was isostatically added on top of the depth-to-basement. The sediment loading approach used here is similar to procedures used by Crough [57] and Sykes [58], was adopted using a multicomponent sediment layer with varying sediment densities calculated from a linear extrapolation of published sediment densities [35]. Once the sediment loading was added, the reconstruction of the open ocean regions was considered complete ( Figure 5), and ready for merging with shelf-slope-rise structures, and adding LIPs, plateaus and seamounts, as will be described next.  The parts of the ocean basins that occupy the transitional zone between oceanic crust and the emergent continental crust are termed here as "shelf-slope-rise" regions. These regions typically extend from the boundary of open ocean regions to the coastline. The three-parameter "shelf-sloperise" model from Goswami [35], based on analysis of various modern shelf-slope-rise structures at active and passive margin regions from ETOPO1, along with their corresponding sediment thicknesses taken from [51], was used to reconstruct the C-T "shelf-slope-rise" regions, extending from the maximum extent ocean crust to the coastlines ( Figure 6).  The parts of the ocean basins that occupy the transitional zone between oceanic crust and the emergent continental crust are termed here as "shelf-slope-rise" regions. These regions typically extend from the boundary of open ocean regions to the coastline. The three-parameter "shelf-sloperise" model from Goswami [35], based on analysis of various modern shelf-slope-rise structures at active and passive margin regions from ETOPO1, along with their corresponding sediment thicknesses taken from [51], was used to reconstruct the C-T "shelf-slope-rise" regions, extending from the maximum extent ocean crust to the coastlines ( Figure 6). The parts of the ocean basins that occupy the transitional zone between oceanic crust and the emergent continental crust are termed here as "shelf-slope-rise" regions. These regions typically extend from the boundary of open ocean regions to the coastline. The three-parameter "shelf-slope-rise" model from Goswami [35], based on analysis of various modern shelf-slope-rise structures at active and passive margin regions from ETOPO1, along with their corresponding sediment thicknesses taken from [51], was used to reconstruct the C-T "shelf-slope-rise" regions, extending from the maximum extent ocean crust to the coastlines ( Figure 6). The paleo-shoreline model used in this reconstruction was obtained from the C-T paleo-DEM obtained from the PALEOMAP Project [41]. In this model, a zero meter coastline was extracted from the paleo-DEM to be used along with the EB08 ocean age grid map. When the extracted continental shorelines from the PALEOMAP Project reconstructions were incorporated into the EB08 CT age grid, there was a very small difference in the rotation of the continental South America block. This was manually adjusted to match with the EB08 grid.
The construction of the shelf-slope model is a combination of geographic information system and hand drawing activity. The process starts with drawing perpendicular lines to the paleocoastline. Once the normal distance between the coastline and the shoreline has been defined, the nearest ocean crust is determined, and the shelf, slope and rise [35] dimensions are marked.
Accordingly, as detailed in [35] the C-T open-ocean regions ( Figure 5) are then merged with C-T shelf-slope-rise regions ( Figure 6). To accomplish this merging, map-based operations such as computing distances between locations were carried out using ArcGIS, and local calculations with interpolation were carried out in Matlab. The merged C-T base reconstructed bathymetry is shown in Figure 7.  The paleo-shoreline model used in this reconstruction was obtained from the C-T paleo-DEM obtained from the PALEOMAP Project [41]. In this model, a zero meter coastline was extracted from the paleo-DEM to be used along with the EB08 ocean age grid map. When the extracted continental shorelines from the PALEOMAP Project reconstructions were incorporated into the EB08 CT age grid, there was a very small difference in the rotation of the continental South America block. This was manually adjusted to match with the EB08 grid.
The construction of the shelf-slope model is a combination of geographic information system and hand drawing activity. The process starts with drawing perpendicular lines to the paleo-coastline. Once the normal distance between the coastline and the shoreline has been defined, the nearest ocean crust is determined, and the shelf, slope and rise [35] dimensions are marked.
Accordingly, as detailed in [35] the C-T open-ocean regions ( Figure 5) are then merged with C-T shelf-slope-rise regions ( Figure 6). To accomplish this merging, map-based operations such as computing distances between locations were carried out using ArcGIS, and local calculations with interpolation were carried out in Matlab. The merged C-T base reconstructed bathymetry is shown in Figure 7. The paleo-shoreline model used in this reconstruction was obtained from the C-T paleo-DEM obtained from the PALEOMAP Project [41]. In this model, a zero meter coastline was extracted from the paleo-DEM to be used along with the EB08 ocean age grid map. When the extracted continental shorelines from the PALEOMAP Project reconstructions were incorporated into the EB08 CT age grid, there was a very small difference in the rotation of the continental South America block. This was manually adjusted to match with the EB08 grid.
The construction of the shelf-slope model is a combination of geographic information system and hand drawing activity. The process starts with drawing perpendicular lines to the paleocoastline. Once the normal distance between the coastline and the shoreline has been defined, the nearest ocean crust is determined, and the shelf, slope and rise [35] dimensions are marked.
Accordingly, as detailed in [35] the C-T open-ocean regions ( Figure 5) are then merged with C-T shelf-slope-rise regions ( Figure 6). To accomplish this merging, map-based operations such as computing distances between locations were carried out using ArcGIS, and local calculations with interpolation were carried out in Matlab. The merged C-T base reconstructed bathymetry is shown in Figure 7.  This reconstructed C-T bathymetry was then updated with known oceanic LIPs and seamounts (Figures 2b and 3b), for which GPlates (www.gplates.org) was used to rotate the modern distribution of the LIPs and seamounts back to C-T time. This updated reconstructed bathymetry was then smoothed and cleaned in Matlab using kriging interpolation. The final map of C-T reconstructed bathymetry is presented in Figure 8.
This reconstructed C-T bathymetry was then updated with known oceanic LIPs and seamounts (Figures 2b and 3b), for which GPlates (www.gplates.org) was used to rotate the modern distribution of the LIPs and seamounts back to C-T time. This updated reconstructed bathymetry was then smoothed and cleaned in Matlab using kriging interpolation. The final map of C-T reconstructed bathymetry is presented in Figure 8.
In sum, the final products are global maps of C-T depth-to-basement at 0.1° × 0.1° resolution (Figure 4), C-T paleobathymetry for the open oceans with an isostatically adjusted multicomponent sediment layer at 0.1° × 0.1° resolution ( Figure 5), C-T paleobathymetry with reconstructed continental shelf-slope-rise structures ( Figure 7) and a realistic combined C-T paleobathymetry with known LIPs and seamounts, continental shelf-slope-rise structures, and isostatically adjusted multicomponent sediment layer at 0.1° × 0.1° resolution (Figure 8).

Results
Here, the general nature of the C-T reconstructed bathymetry is discussed. One important point is that this bathymetric reconstruction was updated only with known LIPs and seamounts. As is clear from Figures 2b and 3b, only certain areas in the central Pacific, Indian, proto-Caribbean and the North Atlantic Ocean therefore could be updated with seamounts and LIPs.

Reconstructed C-T Shelf-Slope-Rise Regions
The reconstructed C-T bathymetry ( Figure 8) indicates that the global distribution of the continental blocks was more clustered than the modern world. North America and Eurasia were still connected, a young north Atlantic was spreading, and the central and southern Atlantic was in a juvenile state. In the south, Indian Ocean was also in its infancy while in between Indian peninsula and Eurasia was huge and mature Tethys Ocean. This unique paleogeography was quite different from the modern world and consequently influenced the C-T shelf-slope-rise regions. Table 2 presents a comparison between the percentages of the global area for continents, oceans, continental shelf, slope and rise. Of the first order, during the C-T time, the area of the land was about 4% less of the global total area than the present world. Here, land refers to the part of the continental plate that is not submerged. Comparing with modern land areas, it is almost 14% less. All of the continental blocks were separated at this point of time, as the rifting of the supercontinent Pangea  Figure 7) and a realistic combined C-T paleobathymetry with known LIPs and seamounts, continental shelf-slope-rise structures, and isostatically adjusted multicomponent sediment layer at 0.1 • × 0.1 • resolution (Figure 8).

Results
Here, the general nature of the C-T reconstructed bathymetry is discussed. One important point is that this bathymetric reconstruction was updated only with known LIPs and seamounts. As is clear from Figures 2b and 3b, only certain areas in the central Pacific, Indian, proto-Caribbean and the North Atlantic Ocean therefore could be updated with seamounts and LIPs.

Reconstructed C-T Shelf-Slope-Rise Regions
The reconstructed C-T bathymetry (Figure 8) indicates that the global distribution of the continental blocks was more clustered than the modern world. North America and Eurasia were still connected, a young north Atlantic was spreading, and the central and southern Atlantic was in a juvenile state. In the south, Indian Ocean was also in its infancy while in between Indian peninsula and Eurasia was huge and mature Tethys Ocean. This unique paleogeography was quite different from the modern world and consequently influenced the C-T shelf-slope-rise regions. Table 2 presents a comparison between the percentages of the global area for continents, oceans, continental shelf, slope and rise. Of the first order, during the C-T time, the area of the land was about 4% less of the global total area than the present world. Here, land refers to the part of the continental plate that is not submerged. Comparing with modern land areas, it is almost 14% less.
All of the continental blocks were separated at this point of time, as the rifting of the supercontinent Pangea continued, though they were not far from each other. Accordingly, the ocean area was 4% larger in terms of total global area, and 5.66% greater the present-day world. These differences may be attributed to the growth in the size of the continents due to collisional tectonics [59].
Oceanic eustasy coupled with terrigenous sediment supply from newly separated continental interiors contributed to the development of mature shelf-slope-rise regions. As discussed in [35], typically the extent of these mature shelf-slope-rise regions extends beyond the maximum extent of ocean crust towards the shoreline. During C-T time, besides the major continental blocks including Africa, America, Eurasia, Antarctica, and India separating from the main Pangean aggregated supercontinent, many small-scale tectonic events took place that contributed to further expansion of shelf-slope-rise regions. For example, Madagascar separated from India, and the Kohistan-Ladakh Arc was rapidly approaching the northern boundary of the Indian peninsula [60]. In the northern Tethys between the North Atlantic and NW Tethys there were several smaller continental blocks. Surrounding both big and small continental blocks, development of shelf-slope rise structures was critical. The presence of an equatorial-tropical seaway connecting the Tethys to the Atlantic oceans, and the shallow but significantly wide Western Interior Seaway connecting the Gulf of Mexico to the Arctic Ocean, and a shallow African Seaway all added to the larger fraction of the shelf-slope-rise area during C-T time. Table 2. Comparison of area between modern and C-T oceans. As summarized in Table 2, the total area occupied by all of the shelf-slope-rise regions was about 22.3% of the global area. Just shelf alone was 8.60% of the global area. Compared to modern combined shelf-slope-rise area, C-T shelf-slope-rise area was larger by more than 30%. The C-T shallow shelf was 8.60% of global area compared to 5.68% for the modern world. In other words, this was more than 60% of modern shelf area. Considering the continental slope and rise together the trend is not so similar, and was about 16% larger during C-T time. In sum, the most significant difference between modern and C-T oceans was in the area of the shelf region. Consequently, from a bathymetric point of view, this additional shallow shelf area and intermediate depth slope-rise area gave the C-T oceanic world a unique dimension.

Reconstructed C-T Open Ocean Region
Beyond the shelf-slope-rise region is the open ocean region. According to Table 2 Positive bathymetric values were changed to 0 using kriging type interpolation methods. Figure 9 shows the distribution of such points on top of the Interpolated 94 Ma reconstructed bathymetry (Figure 8). In this regard, one of the final processing steps involved kriging interpolation to eliminate outlier calculation points. The total number of such points (19,285) was very small compared to the total number of points that define the entire C-T open ocean (4,498,703), i.e., less than 0.29%. Positive bathymetric values were changed to 0 using kriging type interpolation methods. Figure 9 shows the distribution of such points on top of the Interpolated 94 Ma reconstructed bathymetry (Figure 8).  Table 3 shows the statistics for modern and C-T ocean floor age, and the various steps of reconstructed C-T bathymetry in contrast to modern bathymetry (obtained from ETOPO1). The mean age of C-T ocean crust was about 40 million years (Myr) whereas the modern ocean crust has a mean age of about 63 Myr. Thus, on average, C-T oceanic crust was younger than the modern equivalent. The median ages corroborate the same view with C-T being about 36 Myr while the modern ocean is about 55 Myr. Looking at the mode value for C-T oceanic crust a drastic difference is revealed. While the age of modern oceanic crust has most common value of 0 Myr, for C-T oceanic crust it was 190 Myr. Thus C-T oceanic crust had overall older crustal age than the modern ocean. The overall statistics for the C-T reconstructed bathymetry are also provided in Table 3 and Supplementary Figure S3. The output of the plate model, listed as C-T depth-to-basement has maximum and minimum depths of 5555 m and 2663 m, respectively; the average depth is 4465 m. With loading from modeled sediment layers, these values change to 4920, 1385, and 4025 m, respectively. When the shelf-slope-rise wedge model is added to the open ocean reconstruction to obtain the C-T base bathymetry the signature changes to 5760 m, 677 m (above sea level) and 3305 m respectively for maximum, minimum and average depths (Supplementary Figure S4).  Table 3 shows the statistics for modern and C-T ocean floor age, and the various steps of reconstructed C-T bathymetry in contrast to modern bathymetry (obtained from ETOPO1). The mean age of C-T ocean crust was about 40 million years (Myr) whereas the modern ocean crust has a mean age of about 63 Myr. Thus, on average, C-T oceanic crust was younger than the modern equivalent. The median ages corroborate the same view with C-T being about 36 Myr while the modern ocean is about 55 Myr. Looking at the mode value for C-T oceanic crust a drastic difference is revealed. While the age of modern oceanic crust has most common value of 0 Myr, for C-T oceanic crust it was 190 Myr. Thus C-T oceanic crust had overall older crustal age than the modern ocean. The overall statistics for the C-T reconstructed bathymetry are also provided in Table 3 and Supplementary Figure S3. The output of the plate model, listed as C-T depth-to-basement has maximum and minimum depths of 5555 m and 2663 m, respectively; the average depth is 4465 m. With loading from modeled sediment layers, these values change to 4920, 1385, and 4025 m, respectively. When the shelf-slope-rise wedge model is added to the open ocean reconstruction to obtain the C-T base bathymetry the signature changes to 5760 m, 677 m (above sea level) and 3305 m respectively for maximum, minimum and average depths (Supplementary Figure S4). The minimum depth of 677 m above sea level is the result of outliers generated by the mathematical calculations. These outliers were later corrected using kriging type interpolation. The number of these outliers were very small compared to the scope of the entire dataset, (19,285 out of 4,498,703, which is less than 0.29%). The addition of oceanic LIPs and seamounts (both before and after interpolation) did not change the result of fixing the values programmatically for these 19,285 points. These erroneous points (Supplemental Figure S1) are mostly clustered, distributed in low to high latitudes of the Northern Hemisphere, and cover 0.031% of global C-T area (0.04% of global C-T ocean area).
Modern ocean bathymetric statistics derived from ETOPO1 are comparable to the C-T statistics. Although the modern ocean has much deeper trenches, the C-T ocean has a maximum depth of 5758 m, with an average depth of 3233 m, which is surprising similar to the modern (95.44%), which is 3388 m. The median depths of the modern and C-T oceans are also significantly similar, 3863 and 4005 m (103.7%), respectively. The difference may be attributed to the much older overall crustal age of the C-T ocean (see above). The most striking similarity is that the standard deviation of the OES C-T dataset is 100.88% of modern dataset. Another surprising similarity is in variance of these two datasets is within 2% range (Table 3).

Discussion
The reconstructed global C-T bathymetry presented here is the first of its kind based on the observations from modern ocean bathymetry. The main goal of OES is to reconstruct global paleobathymetry, provided age of the oceanic crust, onto which a generalized shelf-slope-rise structure model (based on present-day geometry) is anchored. The EB08 dataset was adopted for oceanic crustal ages. Along with the age of the oceanic crust for the last 140 Myr, EB08 also provides a version of depth-to-basement calculation and thereby a paleobathymetry (global in nature but without shelf-slope-rise) for each of the last 140 Myr in one Myr increments. The comparison of different statistical parameters between models are shown in Supplementary Figure S5.
Accordingly, the discussion starts with the comparison between EB08 (Supplementary Figure S1) and OES. The depth-to-basement reconstructions of OES and EB08 are similar in terms of minimum, mean, median, and mode depths. The difference is observed in the standard deviation, as the EB08 depth-to-basement has a maximum value which is positive. This is due to the fact that EB08 includes LIPs in the depth-to-basement whereas OES reconstruction does not. In Table 3, C-T OES final interpolated reconstructed bathymetry and its various components are compared with EB08, and SW07 bathymetries, which are very different. While OES bathymetry has a mean depth of 3233 meters, EB08 bathymetry has a mean depth of 4269 m. While compared to OES bathymetry trimmed to EB08 extent the mean value for OES bathymetry is 4070 meters. So overall, EB08 reconstructed bathymetry is on an average 200 m deeper in the open ocean regions. The OES full extent bathymetry is but much shallower than EB08 because of the shelf-slope-rise regions. The major difference in the EB08 bathymetry and the OES bathymetry EB08 extent is clear in variance, the measure of how far each value in the data set is from the mean. The variance for the OES bathymetry dataset (full bathymetry) is extraordinarily analogous to the variance of modern world bathymetry. The OES bathymetry dataset (EB08 extent), compared to EB08, has much smaller variance signifying gradual changes in bathymetry than that in EB08. This is confirmed by the modal and median values which shows that OES bathymetry has shallower median values as well as very shallow depths as the most common values. In EB08, the mode is very high at 5150 m as opposed to 6.5 m in OES. When compared to the OES bathymetry trimmed to EB08 extent, high modal value 4878 m (though much less than EB08, 5150 m) is observed. This exactly spells out the difference between the two bathymetry dataset and highlights the importance of reconstructing the shelf-slope-wedge. This has special significance especially for the C-T time, as for decades the paleoclimate community has argued for extensive shallow epicontinental seas [36,[62][63][64][65][66][67][68][69][70][71][72][73][74][75][76] for example, to explain oceanic anoxic events. The OES C-T Ocean compared to Modern Ocean, has 4.75% more area (ocean area) less than 1000 m deep (maximum extent of combined euphotic and disphotic zone depth). For areas shallower than 500 m and 200 m, the OES C-T Ocean has 3.69% and 1.21% more compared to the Modern Ocean (Supplementary Table S1). The corresponding maps showing the extent of the reconstructed OES C-T ocean with these depths are shown in Supplementary Figure S6. These increased areas of shallower depths (compared to total ocean areas) are encouraging because OES reconstructed bathymetry at least in the case of C-T time, reconstructed larger areas of shallower depths in agreement with previously postulated arguments.
The difference between the depth-to-basement and final bathymetry also sheds some light regarding how much sediment was used in these bathymetric reconstructions, although caution is needed as direct measurement of sediment amount was not involved. Isostatic loading of sediment is complicated owing to density and thickness of sediment layers. The OES model shows a significant difference in the depth-to-basement to final bathymetric reconstruction ( Table 3, OES Difference). Addition of sediment, modeled shelf-slope-rise structures and known seamounts and oceanic LIPs significantly change the final bathymetry with a much larger standard deviation. This is also seconded by the median values of the two models: OES median value is 48 m and EB08 median value is 205 m. This implies that the changes in depth-to-basement to final bathymetry in the OES reconstruction does not involve many abrupt changes (or less abrupt changes) in depth, and rather, follows a natural gradient, which may be more realistic. A high median value, on the other hand, implies that in EB08, there are unexpected abrupt variations in the bathymetric values. The most common value (mode) in the OES difference dataset is 675 m, while the mean is 385 m. In EB08, these values are considerably lower, 490 m and 250 m, respectively.
When compared with the modern ocean bathymetry the OES C-T reconstructed bathymetry reflects interesting similarity. Though it is created using the very different ocean crustal age (see Table 3) than the modern one, the average and median depth, mode, standard deviation and the variance of the modern ocean and C-T ocean are highly comparable (Table 3). Surprisingly, these two datasets (modern and OES C-T) have less agreement with the EB08 reconstructed bathymetry.
Other than statistical comparison, a difference map between the EB08 and OES C-T reconstructions is given in Figure 10. The red tone signifies that EB08 reconstruction is deeper while blue tone signifies the OES reconstruction is deeper. Statistics for this difference map are provided in Table 3, but also computed EB08 extent (mainly open oceans regions) as EB08 does not have shelf-slope-rise structures. Both maximum and minimum values for the difference plots are very large. The maximum value of 5640 m are in the regions where there are underwater LIPs. In EB08, these areas have positive values signifying they are above the water and sometimes even up to 2000+ m. In the OES reconstruction, all known LIPs are added, but they all have an average depth of at least 1500 m. Similarly, the high minimum number is near the outer (oldest) edge of the oceanic crust, where in the OES reconstruction, the shelf-slope-rise model is added to estimate natural coastal landforms.
Much of the open ocean regions in the two models are in agreement and their mean difference is 190 m with a high standard deviation.
Another interesting comparison can be made with the C-T bathymetric model of Sewall [2] (SW07). The SW07 model is available in a 64 × 128 grid size, and here is resized to 1801 × 3601 (Supplementary Figure S2) in order to compare with the OES model. The SW07 model also contains topography; only values less than or equal to zero were considered for comparison ( Figure 11). The statistical comparison between the OES and SW07 models indicates some fascinating similarity. The minimum, mean and median depths are within 10% of each other. The modal depth, however, is very different, in OES bathymetry being very shallow (6.5 m), while that in SW07 is very deep (5040 m). The standard deviations of these two models are 10% (OES) and 20% (SW07). In the SW07-OES difference map (Figure 11 In addition to the Cretaceous paleobathymetric reconstruction presented here, there are other important related developments that are very encouraging. Hochmuth and Gohl [33] are working on improving geodynamic and paleoceanographic aspects related to ocean crustal age with a highresolution plate kinematic model. The alternative seafloor age grid for the South Atlantic Ocean based on a recent high-resolution plate kinematic model [43] is of particular interest due to its connection to the Cretaceous Normal Superchron. Enhanced studies like this will open the door for better understanding of long-offset fracture zones, currently mostly overlooked in paleobathymetric reconstructions.
The C-T transition occurred during the wettest, hottest climates of the Cretaceous greenhouse with assists from greenhouse gases emitted to the atmosphere by unusually large mid-ocean ridge systems and mantle plume magmatism [77]. The opening of the critical N-S Equatorial Atlantic Seaway coincided with a change from an "eddy ocean" to a "thermohaline ocean" [78]. Evidence for this change is provided by North Atlantic nannofossils which evolved from homogeneous Cenomanian distributions with cosmopolitan ranges to heterogeneous Turonian distributions with restricted ranges [40].
The Earth's equator-to-pole temperature gradient is a fundamental climatological property related to equator-to-pole heat transport efficiency, ocean circulation, hydrological cycling and atmospheric chemistry [79]. The meriodional thermal gradient may have been at its lowest in the Cretaceous during the C-T transition. Unfortunately, Cretaceous paleoclimate modeling with GCMs has consistently underestimated the warm polar temperatures that are indicated by fossil and geochemical data (e.g., [74,80,81]).
Realistic paleobathymetries are among the improvements that can be tested in GCM paleoclimate simulations. The goal of this project has been to address this by producing a more realistic ocean floor with continental shelf-slope-rise, underwater ridges, trenches, seamounts, LIPs and a corresponding ocean bottom roughness. In particular, seafloor roughness and topography exert critical controls on vertical ocean mixing and circulation [82,83]. Will our paleobathymetry provide clues for understanding the eddy-dominated Cenomanian ocean? What would cause ocean reorganization in the Turonian-development of a less rough seafloor due to more rapid seafloor spreading [42,83].

Conclusions
The ongoing goal of OESbathy reconstruction is to provide a means for reconstructing realistic paleobathymetries for geologic times that can be tied directly to ocean crustal age. The innovation of OESbathy is the generalized shelf-slope-rise model to provide bathymetric interpolation across the transitional zone between open ocean and coastlines; this model is described in detail for the test case modern day OESbathy 1.0 reconstruction [35]. Here OESbathy has been extended to Cenomanian-Turonian (C-T) boundary time (94 Ma), and updated with known oceanic LIPs and seamounts. The main results are as follows: Figure 11. Difference map between the SW07 and OES C-T reconstructions.
In addition to the Cretaceous paleobathymetric reconstruction presented here, there are other important related developments that are very encouraging. Hochmuth and Gohl [33] are working on improving geodynamic and paleoceanographic aspects related to ocean crustal age with a high-resolution plate kinematic model. The alternative seafloor age grid for the South Atlantic Ocean based on a recent high-resolution plate kinematic model [43] is of particular interest due to its connection to the Cretaceous Normal Superchron. Enhanced studies like this will open the door for better understanding of long-offset fracture zones, currently mostly overlooked in paleobathymetric reconstructions.
The C-T transition occurred during the wettest, hottest climates of the Cretaceous greenhouse with assists from greenhouse gases emitted to the atmosphere by unusually large mid-ocean ridge systems and mantle plume magmatism [77]. The opening of the critical N-S Equatorial Atlantic Seaway coincided with a change from an "eddy ocean" to a "thermohaline ocean" [78]. Evidence for this change is provided by North Atlantic nannofossils which evolved from homogeneous Cenomanian distributions with cosmopolitan ranges to heterogeneous Turonian distributions with restricted ranges [40].
The Earth's equator-to-pole temperature gradient is a fundamental climatological property related to equator-to-pole heat transport efficiency, ocean circulation, hydrological cycling and atmospheric chemistry [79]. The meriodional thermal gradient may have been at its lowest in the Cretaceous during the C-T transition. Unfortunately, Cretaceous paleoclimate modeling with GCMs has consistently underestimated the warm polar temperatures that are indicated by fossil and geochemical data (e.g., [74,80,81]).
Realistic paleobathymetries are among the improvements that can be tested in GCM paleoclimate simulations. The goal of this project has been to address this by producing a more realistic ocean floor with continental shelf-slope-rise, underwater ridges, trenches, seamounts, LIPs and a corresponding ocean bottom roughness. In particular, seafloor roughness and topography exert critical controls on vertical ocean mixing and circulation [82,83]. Will our paleobathymetry provide clues for understanding the eddy-dominated Cenomanian ocean? What would cause ocean reorganization in the Turonian-development of a less rough seafloor due to more rapid seafloor spreading [42,83].

Conclusions
The ongoing goal of OESbathy reconstruction is to provide a means for reconstructing realistic paleobathymetries for geologic times that can be tied directly to ocean crustal age. The innovation of OESbathy is the generalized shelf-slope-rise model to provide bathymetric interpolation across the transitional zone between open ocean and coastlines; this model is described in detail for the test case modern day OESbathy 1.0 reconstruction [35]. Here OESbathy has been extended to Cenomanian-Turonian (C-T) boundary time (94 Ma), and updated with known oceanic LIPs and seamounts. The main results are as follows: • Parameterized shelf-slope-rise structures result in a much shallower mean ocean depth than other available C-T bathymetry reconstructions. • C-T ocean crust was on average 40.29 million years old whereas modern ocean crust has a mean age of 63.39 million years.

•
The OES C-T ocean has an average depth of 3262 m, while the modern ocean average depth is 3388 m. Despite that the OES C-T Ocean is about 125 m shallower, the overall ocean bathymetric pattern matches with the modern bathymetric values as both have similar standard deviations (within 1%) despite significantly different maximum and average depths. • Seafloor roughness is provided from mid-ocean ridges, and oceanic LIPS and sea mounts rotated to their C-T positions (using known features).
The following improvements are being planned: (1) The OES sediment model will be validated and corrected using data from DSDP/ODP/IODP. (2) The OES shelf-slope-rise model will be revisited with a more detailed parameterization of modern active and passive margins. (3) We are working on the reconstruction of ocean bottom surface fractures, especially those near the mid-ocean ridge areas, and experimenting with statistical models to simulate seafloor roughness. (4) An OES reconstruction of Paleocene-Eocene Thermal Maximum (55 Myr) paleobathymetry is now underway.
Supplementary Materials: The following are available online at www.mdpi.com/2076-3263/8/1/21/s1; Figure S1: EB08 C-T bathymetry; Figure S2: SW07 [2] C-T bathymetry (from the total land ocean topography bathymetry, only the values zero or below have been plotted here); Figure S3: Comparative statistics of modern and C-T ages of the ocean crust (in Ma). The minimum is not shown as both have zero as minimum ocean crust age; Figure S4: Comparison of the three models. OES: Interpolated OES CT Bathymetry, OES*: Interpolated OES CT Bathymetry trimmed to EB08 extent, EB08: EB08 CT Bathymetry, and SW07: SW07 CT Bathymetry; Figure S5: Comparison between different of OES C-T bathymetry development.
Step 4: OES Bathymetry with known LIPs and Seamounts, and Step 5: OES Interpolated Bathymetry; Figure Table S1: Comparison of shallow ocean areas between Modern and C-T oceans.