A generalised approach to reconstructing geographical boundary conditions for palaeoclimate modelling

Introduction Conclusions References


Introduction
As the interest grows in assessing how the climate may change under increasing anthropogenic forcing, the field of climate modelling is steadily growing and improving.These models are mainly tested by making hindcasts of the present day and pre-industrial climate.Anthropogenic forcing may lead to future climate states that are beyond the range for which models can be validated with observations.Reconstructing the palaeoclimate may help to increase the confidence in predictions by simulating episodes in Earth's history with climatic conditions that were drastically different to that of today.While the number of measurements is generally much more limited than for the present day, a considerable amount of proxies is available to constrain the conditions of several such periods, like the Paleocene-Eocene thermal maximum, Eocene-Oligocene transition and Eemian.The equable climate problem (Huber and Caballero, 2011) is one good example showing that models still lack some essential dynamics to explain how the meridional atmospheric (or sea surface; Bijl et al. (2009)) temperature gradient can become so flat in a greenhouse climate.It is therefore essential to try and simulate the past climate with a similar amount of detail as is done for the present to better understand the climate system when it is forced far beyond its present state.
Simulating deep time palaeoclimate is challenging, since the applied boundary conditions can be poorly constrained.Between reconstructions, differences in e.g.lat/lon coordinates, heights, coastlines and basin depths are are often considerable.Next to the atmospheric composition, solar insolation and land surface properties, the global geography can be of great importance (Barron and Washington, 1984).The breakup and collision of continents likely drove major climate transitions in the past.For example, the India-Eurasia collision caused a large scale uplift that resulted in regional climate changes but also affected the global circulation and hydrological cycle (Ramstein et al., 1997;Dupont-Nivet et al., 2007, 2008;Caves et al., 2014).Furthermore, the opening of ocean gateways such as the Tasmanian Gateway impacted ocean circulation and thus the redistribution of heat, resulting in both changed surface and deep ocean temperatures (Bijl et al., 2013;Sijp et al., 2014).Many of these themes are the subject of ongoing studies and continuous acquisition of new data, either of palaeo-altimetry/bathymetry or from tectonic reconstructions.It is therefore crucial that new relevant insights can be incorporated easily and swiftly into boundary conditions used for palaeoclimate modelling.
Implementing different geographies into climate models is not without issues; the increasing but limited model resolution obscures certain potentially critical features (such as narrow ridges, mountain ranges and straits).Additionally, plate-tectonic reconstructions have their own uncertainties and are continuously being expanded and improved.This causes large discrepancies between model simulations because different studies use reconstructions of varying age and approach.A simple but striking example is given by Gasson et al. (2014), who used the output of several coupled global atmosphere-ocean models to force an ice sheet model simulating the Antarctic glaciation during the Eocene-Oligocene transition (∼34 Ma).Their main conclusion states that it is nearly impossible to determine the CO 2 threshold for the inception of ice sheets since the models use such different topographies.Geographical reconstructions used in palaeoclimate modelling are often based on older reconstructions (such as Markwick (1998); Sewall et al. (2000)) that are adjusted for a specific case.
These reconstructions may thus be based partially on information that is outdated and are usually only made for a unique set of time frames, which greatly limits their application in other studies Since small-scale features can be crucial in the climate system, they always need to be considered and may require additional manual changes.Finally, developments in coding and computer hardware allow for model simulations with a significantly higher resolution that in turn require more detailed boundary conditions.This calls for a new approach where the latest geological reconstructions can directly be incorporated together with improving constraints on palaeoaltimetry (e.g.Dupont-Nivet et al. (2008)) and translated into a gridded topography and bathymetry map.In this paper, a method is presented that allows for the creation of geographical reconstructions based on the most recent geological insights with relative ease.Since there is a relative scarcity in middle and late Eocene reconstructions, the composition of a 38Ma palaeogeography will be presented and serve as an example of the general workflow.The set of files (and downloads) required to reproduce this procedure for 38Ma, as well as other time frames, is provided in the supplementary material along with a Readme for guidelines.
Although somewhat similar workflows have been presented in the past (Markwick and Valdes, 2004;Herold et al., 2008), the innovative aspects here are the approach of directly incorporating a platetectonic reconstruction, the freedom of choosing a different plate-tectonic model framework and the ability to consider multiple time frames.

Plate tectonic reconstruction
The backbone of our method is a plate-tectonic model that provides a first-order, 200Ma to present day reconstruction of palaeogography (Seton et al., 2012).This model is provided in the free plate reconstruction software GPlates (www.gplates.org,Boyden et al. (2011)), which is widely used in tectonic research.The reconstruction of Seton et al. (2012) serves as a basis, and can be updated when new reconstructions become available.Because we aim to study palaeoclimate, we place our reconstruction in a palaeomagnetic reference frame (Torsvik et al., 2012;van Hinsbergen et al., 2015).Another reference frame can be used by simply changing the rotation file of the plate-tectonic model, only the ocean depth grid needs to be consistent with the latter.Starting from the presentday topography, the outlines of land, continental shelves and shallower regions such as intra-oceanic ridges and plateaus are defined.Using the reconstruction model, these features are then shifted towards their positions at a chosen time frame (here: 38Ma) and exported as Shapefiles (.shp).The latter are then transformed into numerical grids and stored as a NetCDF file using a simple Matlab script for easy and unambiguous interpretation by most computational programs.An example showing the gridded masks for the middle-late Eocene is shown in Figure 1.
For its application as a climate model boundary condition, the final geographical reconstruction should be stored as a gridded (latitude-longitude) topography-bathymetry dataset.To start with, ocean bathymetry is estimated using the present day ocean floor age grid provided by Müller et al. (2008a) in combination with the plate tectonic model.Using the latter, the agegrid is shifted into its position at the considered time frame (38Ma here) to be compatible with the rest of the reconstruction.The position and age (corrected for the plate motions that have occurred since then) of oceanic plates in the model is then translated into a bathymetry map using the general depth-age relationship for oceanic crust.Here, we use an average of the relations given for the main ocean basins (Atlantic, Pacific and Indian Ocean) treated by Crosby et al. (2006): where Z is the basin depth in m and A the ocean floor age in millions of years.
The method above can only be applied to regions where the ocean floor is still preserved today, leaving large gaps where oceanic crust has subducted since the reconstruction time (e.g. the western and southeastern Pacific Ocean).An estimated palaeo-age grid was already translated into a rastered geo-referenced depth to basement dataset by Müller et al. (2008a).The latter also provides a reconstructed bathymetry, including sediments, but this makes it incompatible with other areas where only the age-depth relationship is used.Furthermore, sediments are only significant in a few locations where the ocean floor is still very deep (generally below 5000m), resulting usually only in small changes.This is even more so, when the bathymetry grid is translated in depth layer masks that are used in numerical models, as they rarely go beyond 5-5,5km.One exception to this is the Neotethyan Ocean where up to 1500m of sediment is added in Müller et al. (2008a), but this region is treated specifically here resulting in a similar final depth (4,5km).Finally, the age-depth relationship used here (Crosby et al., 2006) is very similar to the one used to reconstruct the depth to basement grids (Stein and Stein, 1992) as differences are generally around 50m or less.
The Müller et al. (2008a) grids are provided in a hotspot reference frame and thus need to be converted to the palaeomagnetic reference frame used here in the scope of climate modelling (Moreau et al., 2007;van Hinsbergen et al., 2015).We applied this conversion by fixing the entire bathymetry grid to Africa which is the central pole in both of the reconstruction frameworks.First, the movement of Africa in the past 38Ma is reversed by using an adapted version of the Seton et al. (2012) rotation file such that it moves back to its present position starting from where it was at 38Ma instead of the opposite.The bathymetry map (defined at 38Ma) is then shifted to the configuration it would have at the present day, disregarding the individual plate movements relative to Africa.Finally, the same is done using the (unaltered) rotation file of the (Torsvik et al., 2008) palaeomagnetic framework, moving the bathymetry grid back to its configuration at 38Ma as if it were made using this reconstruction model.After doing so, the adapted Müller et al. (2008a) depth grid fits within the continental margins of the plate-tectonic reconstruction.Differences between the rotational frameworks used may lead to discrepancies in the reconstruction, but since oceanic plateaus and ridges are treated separately in our approach this does not result in any problems over at least the Cenozoic.For earlier periods, more specific adjustments may be necessary to cope with possible inconsistencies.In addition to the deep ocean, shallow basins and plateaus are adjusted to their present-day depth at the is taken here, values between 50 and 500m are found in literature) as they mark the outer margins of continental crust, that are assumed to be fairly robust over time.It is important to consider here that the shape and depth of continental shelves may be subject to change, depending on the considered time frame (Sømme et al., 2009).Since the middle-late Eocene had a greenhouse climate, a flat and relatively shallow shelf is implemented here but this can easily be done differently for other periods.
The next, more challenging step is creating a realistic land topography for a chosen time frame.
Good estimates of altimetry are only available for a limited number of locations and time intervals, making the reconstruction of palaeotopography a difficult task that requires a lot of manual input.
For example, the global Eocene topography by Markwick (1998) is still used and referred to regularly in palaeoclimate modelling.Our approach here is to first consider the present day topography (using the ETOPO1 dataset, http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1)and adjust it to the 38Ma configuration using the plate tectonic model and the land mask shown in Figure 1.Later on, specific topographical features can be corrected manually, depending on the specific demands of the application and the amount of information that is available.Since the topography is coupled to the plates and shifted within GPlates, there are no issues of area preservation or distortion.The resulting first stage of the reconstruction is given in Figure 2.

Generic adjustments
The geography obtained at this stage (Figure 2) is not yet accurate and contains a number of gaps that remain to be filled: mainly collision zones and transition regions between continental and oceanic crust.Especially the land topography needs to be adjusted so that it better matches palaeoaltimetry estimates for the considered time frame (38Ma).Although glaciers possibly existed on Antarctica (Scher et al., 2014) and Greenland (Eldrett et al., 2007), the presence of continental-scale ice sheets in the Eocene is unlikely (Deconto et al., 2008).Therefore, to correct for the effect of isostasy the present-day bedrock topography is raised by 1/3 of the current ice thickness.The changes mainly affect Greenland and Antarctica, although for the time window chosen here a highly detailed palaeotopography of Antarctica was derived by Wilson et al. (2012), considering both isostatic rebound and sediment restoration with respect to the present-day topography.The latter, however, assumes the Antarctic continent to be at its present-day location which does not fit into this reconstruction.Similar to how the topography was treated above, the reconstructed Antarctic palaeotopography is fixed to the plate polygons defined in the plate-tectonic model and adjusted towards its 38Ma configuration.In practice, only the area above sea level is incorporated in the reconstruction and combined with the previously determined bathymetry.This is done to minimise possible discrepancies between regions of the ocean floor that are reconstructed using different methods.
Although most of the ocean floor has been reconstructed, a large gap is still present in the Neotethyan region between Africa, Eurasia, India and Indochina at 38Ma.While in the circum-Pacifc area the crust lost to subduction was entirely oceanic, the former was characterised by mutliple continental collision and subduction episodes (Stampfli and Borel, 2002).The Neotethyan Ocean is introduced manually, based on the incorporation of GPlates reconstructions covering the Mediterranean region by van Hinsbergen et al. (2012van Hinsbergen et al. ( , 2014)), the Arabia-Eurasia collision zone by McQuarrie and van Hinsbergen (2013), the Arabia-India plate boundary zone by Gaina et al. (2015) and the India-Asia collision zone by van Hinsbergen et al. (2011van Hinsbergen et al. ( , 2012)).A deep ocean basin with a depth of 4500m is assumed, which is close to the value used in Müller et al. (2008a), representing the Neotethyan Ocean at 38Ma.
Next, we adjust the elevation of several mountain ranges towards values that are more representative for the considered time frame.Since most of the current elevated regions underwent a considerable part of their uplift after the late Eocene (such as the Andes, Himalayas and East African Rift) they will have to be lowered, while eroded regions may require some uplifting (for instance large river catchments).More specific changes for the example considered here will be discussed in more detail in Section 4, the result of the generic changes is first shown Figure 3.As this is clearly the most tedious part of the procedure, it will consume a large part of the time required to acquire a palaeogeography for another time frame.However, time-specific adjustment masks are created as polygons within GPlates and fixed to the plate movements.Consequently, only their extents and effects will need to be adjusted and some masks may need to be removed or added, which facilitates this step significantly.When available, the incorporation of a digital palaeo-elevation model, such as the one presented by Markwick and Valdes (2004) can replace the present day topography so that no further adjustments are needed.

Specific adjustments for the middle-late Eocene
The example shown thus far is largely representative for the middle-late Eocene (38Ma) geography, but now a number of time-frame specific adjustments are made that have not yet been treated in the previous sections.Several major tectonic events took place prior to and during the Eocene-Oligocene transition that had an impact on the global climate.A brief overview of these events and how their potential impacts are implemented specifically is given below and a map of the masks used to apply the various changes is provided in Figure 4.The combined effect of these adjustments is indicated by the highlighted areas in Figure 5, note that the effects of the land mask shown (in white) in Figure 4 have not yet been implemented.The latter is used to correct the land boundaries after the manipulations discussed in Section 5.
The collision between Greater India and Eurasia is causing profound changes to the surrounding geography (van Hinsbergen et al., 2012).The Neotethyan Ocean narrows due to the combined north-ward propagation of Africa and India, leaving the Mediterranean Sea and the northern Indian Ocean.
Further north, a major intra-continental seaway known as the Paratethys covered much of central Europe and a considerable part of western Asia.This was mainly a shallow basin except for the South Caspian and Black Sea areas, as it was mostly underlain by unextended continental crust.The extent of the Paratethys as well as the distribution of the surrounding continents is determined using the GPlates reconstructions mentioned in Section 3. The Neotethyan subduction-collision systems also create a large-scale uplift in southern Asia with primarily the formation of the Himalaya mountain range and the Tibetan Plateau.To the north and west, this uplift forces the westward retreat of the relatively shallow Paratethys leaving several isolated basins (Meulenkamp and Sissingh, 2003;Bosboom et al., 2014a).In addition, the formation of a rift between Antarctica and South America causes the opening of the Drake Passage with the first formation of oceanic crust at around 35 Ma (Livermore et al., 2005).A similar process is taking place south of Australia with the formation of the Tasmanian Gateway.While its opening already occurs around 50Ma, the Tasmanian Gateway initially remains shallow.Coincident with the transition from a slip-boundary towards North-South spreading between Australia and East Antarctica (Bijl et al., 2013), a profound deepening phase occurs from 35.5Ma onwards.The deepening and widening of the Tasmanian gateway is further complicated by a northward extension of the Antarctic shelf and an elevated Tasman Rise which both sank over time (Stickley et al., 2004;Close et al., 2009).Dependent on which absolute plate reconstruction frame is used, the opening of the Tasmanian Gateway can occur either to the north or south of 60 • S paleolatitude.This can have important consequences for the direction of a possible initial throughflow and again stresses the importance of correct paleogeographies in paleoclimate modelling.The related northward motion of Australia combined with the rotation and rifting taking place around Indochina also causes a gradual blocking of the Indonesian Throughflow during the Oligocene (Hall, 2002).
As was pointed out before, some smaller geographical features have also been treated individually.
Mountain ranges that are believed to have experienced major uplifts since the Eocene are mainly the Andes, Himalayas, Sierra Nevada (USA), East African Rift and Middle East, which are all lowered significantly.The southern Tibetan Plateau was lowered more than the Himalayas since a narrow yet high ridge was already present while the adjacent plateau was still much lower (van Hinsbergen et al., 2012).The northern part of the Tibetan Plateau is changed into lowland that surrounds the retreating Paratethys in the Tarim Basin (Dupont-Nivet et al., 2008;Bosboom et al., 2014b;Sobel et al., 2013).
Related to the changes in Eurasia mentioned above, the narrow Tethys is drawn as a deep ocean while the Paratethys consists mainly of shallow seas with the Turgai Strait almost closed (Bosboom et al., 2014a).The latter forms the connection between the Arctic Ocean and the Paratethys and is an important contributor to the water budget in Central Asia.During the early Eocene, the Arctic Ocean may have been geographically isolated as indicated by freshening reported in Brinkhuis et al. waters (similar to the Indonesian Throughflow today), as it consists of an intra-oceanic volcanic arc since the late Cretaceous (Boschman et al., 2014).The Ninety-East Ridge (Slotnick et al., 2015) was also enhanced (made wider and higher) as it is a pronounced yet narrow feature in the eastern Indian Ocean that would otherwise be obscured by smoothing later on.Regarding the areas discussed above, three masks are applied defining ocean (deep oceanic), sea (deep continental) and shallow basins (shallow continental) implying a depth of 4500m, 1000m and 100m, respectively.

Final modifications and smoothing
All major features of the geological reconstruction are now present but before it can be used as a numerical grid for a climate model simulation, some additional modifications still need to be made.Firstly, to estimate the remaining gaps where the altimetry or bathymetry is still unknown, an interpolation is performed using the inpaintn algorithm developed by Garcia (2010) and Wang et al. (2012).An additional land mask was made (shown in white in Figure 4) containing regions that should be above sea level to avoid creating artificial basins with this technique.The land mask is based on several global reconstructions such as those made by Dr. R. Blakey (https: In the middle-late Eocene, mean sea level was about 100m above that of today (Müller et al., 2008b), so we decrease the elevation (and similarly increase the depth) of all data points by 100m to correct for this.The depths that were implemented before, using the masks shown in Figure 4, were initially corrected so they end up being the actual values stated above after changing the global sea level (such that an area assigned to be e.g.1000m deep does not end up being 1100m instead).Again, the land mask is used to correct the elevation of certain inundated regions back to sea level where land is believed to be present at the considered time.For instance, river basins such as the Amazon, Rio de la Plata, Congo and Mississippi were just carved out deeper after the sea level dropped and other depressions were only created during glacial periods (e.g. the Great Lakes).
Finally, the application of the geographical reconstruction will usually be at a lower resolution than the one used until now (0.1 • ) and the high level of detail in the present-day topography may be misleading with respect to how precise palaeotopography is actually known.The last step is then to apply a squared smoothing mask, performing a 2-dimensional running mean with a 1 • width for the example, of which the result is shown in Figure 6.The mask size can be chosen differently and after smoothing, the global geography data is extrapolated onto the desired output grid (e.g., longitudelatitude) to be stored as a NetCDF file.

Conclusions
We have presented a new method to create global geographical reconstructions that can be used as boundary conditions in palaeoclimate model simulations.Starting from a plate-tectonic reconstruction, the depth of the ocean floor was estimated and combined with a reconfigured present day topography.An overview of the different steps making up the procedure is given in Fig. 7. Reconstructed plate trajectories, a background rotational frame and a set of boundary lines are fed into GPlates to make a global plate-tectonic reconstruction for a particular time frame.Within the same setup, present day topography and ocean floor data are manipulated and shifted to fit the same geographical configuration as the rest of the output.The latter is then used in a MatLab script together with some external input (mostly from literature) to create a global gridded altimetry-bathymetry dataset that is stored as a NetCDF file.This results in a highly detailed dataset that needs some timespecific changes, after which it is smoothed and extrapolated onto the desired computational grid.
Even though we performed a number of manual adjustments, interfacing GPlates with our Mat-Lab script allows for the development of global geographical reconstructions with limited effort.
The primary focus was to use a plate tectonics model not originally designed to create numerical grids as output (but it can use and manipulate them), that can be altered afterwards without having to run through the entire process again (only some of the masks may need to be redefined).Creating a new grid for a different point in time or a new configuration (e.g. using a different rotation framework) should allow the user to make multiple reconstructions in a few days' time.In case of a time slice experiment, one can also create appropriate boundary conditions for each time frame instead of making some minor adjustments to an average geography, to really try and catch the effect of gradual but global changes.The amount of manual handling depends on the demands of a specific application and can consequently be very limited, but the procedure does not inhibit additional adjustments for further detail and correctness.In the example treated here, the output grid is defined as having a 0.1 • resolution but can easily be extrapolated onto a coarser grid, which may require some more elaborate smoothing.
As the procedure starts from a plate tectonic model, most of the features important for modelling  (2012).Herold et al. (2014) have manually added a lot of detail to both the land surface and ocean floor in their reconstruction but here the problem remains that there can be considerable differences in continental positions as the geography is bound to a specific time frame.Details within the contours of an older reconstruction can still be added but changing the relative positions of various land masses is not straightforward.
In summary, the reconstruction method presented here can provide detailed palaeogeographical datasets with relative ease for climate model studies and permits users to customise the result to their specific needs.This should facilitate the communication between geologists and climate modellers in their efforts to improve the quality and consistency of future palaeoclimate simulations.
Figure 7. Flowchart of the procedure; the plate-tectonic reconstruction is made with GPlates, using a set of reconstructed plate trajectories (Seton et al., 2012;van Hinsbergen et al., 2015), a reference rotational frame (Torsvik et al., 2008) and a set of boundary masks delimiting present day coastlines, continental shelf breaks, shallow ocean regions and some regions that need manual adjustments (see Figures 1 and 4).In addition to the reconstruction, present day (PD) topography and bathymetry base files are manipulated using the same rotational frame and adjusted to the configuration at the aimed time frame.The resulting output files as well as some external input are then fed into a MATLAB script which builds up a global geography dataset through several consecutive steps and stores it as a NetCDF file. Acknowledgements.
-This work was carried out under the program of the Netherlands Earth System Science Centre (NESSC) 350

Figure 1 .
Figure 1.Earth reconstruction from the plate tectonic model at 38Ma showing shifted present day land surface (white), continental shelves (blue) and relatively shallow ocean regions (black).Remaining in grey are then the deep ocean and zones that are yet unidentified in the plate-tectonic reconstruction model.

Figure 2 .
Figure 2. First stage of the reconstruction: global geography at 38 Ma using the depth-age relationship, the Müller et al. (2008a) reconstruction, shallow basins, continental shelves and the present-day topography.The areas shaded in light grey are gaps which still have to be addressed.

Figure 3 .
Figure 3. Second stage of the reconstruction: topography and bathymetry with adjustments to polar regions and a lowering of several mountain ranges.Changes can be noted mainly in Greenland, Antarctica, the Neotethyan Ocean, the Andes and the Himalayas.

Figure 4 .
Figure 4. Masks used for manipulation of the reconstruction; white represents land, cyan shallow water, blue sea and black ocean.In southern Asia, green, yellow, orange and red indicate areas that are lowered manually by different ratios.

Figure 5 .
Figure 5. Reconstructed topography and bathymetry with specific adjustments indicated by shaded areas, these include mainly Antarctica, southern Panama, the Amazon and several other river basins, the Hudson Bay and Great Lakes, Neotethyan Ocean and Paratethys Seaway, the Ninety-East Ridge and the Tasman Rise.

Figure 6 .
Figure 6.Final smoothed global topography and bathymetry reconstruction for the middle-late Eocene at 38 Ma.
applications are quite realistic.The width and depth of ocean basins, presence of oceanic ridges and plateaus, position of continents and opening or closing of gateways are all directly incorporated in the final result.The set of boundaries shown in Figure1are coupled to individual plates and can easily be linked to a different reconstruction framework.When comparing the result (Figure6) to other recent reconstructions used in palaeoclimate modelling, several differences quickly become evident.For instance, the global geography used inLadant et al. (2014) is coarser and contains less detail in both land topography and ocean bathymetry.The Southern Ocean gateways are also much wider, which is probably a result of both the large grid spacing used for the climate model and the difference in time frame considered as this reconstruction is based on an Oligocene geographyLefebvre et al.