Geoscientific Model Development A dynamic continental runoff routing model applied to the last Northern Hemisphere deglaciation

We describe and evaluate a dynamical continental runoff routing model for the Northern Hemisphere that calculates the runoff pathways in response to topographic modifications due to changes in ice thickness and isostatic adjustment. The algorithm is based on the steepest gradient method and takes as simplifying assumption that depressions are filled at all times and water drains through the lowest outlet points. It also considers changes in water storage and lake drainage in post-processing mode that become important in the presence of large ice dammed proglacial lakes. Although applicable to other scenarios as well, the model was conceived to study the routing of freshwater fluxes during the last Northern Hemisphere deglaciation. For that specific application we simulated the Northern Hemisphere ice sheets with an existing 3-D thermomechanical ice sheet model, which calculates changes in topography due to changes in ice cover and isostatic adjustment, as well as the evolution of freshwater fluxes resulting from surface ablation, iceberg calving and basal melt. The continental runoff model takes this input, calculates the drainage pathways and routes the freshwater fluxes to the surface grid points of an existing ocean model. This results in a chronology of temporally and spatially varying freshwater fluxes from the Last Glacial Maximum to the present day. We analyse the dependence of the runoff routing to grid resolution and parameters of the isostatic adjustment module of the ice sheet model.


Introduction
The problem of finding the runoff pathways and drainage basins for a given topographic setting is crucial for questions in a wide range of scientific fields ranging from hydrology to climatology and environmental studies, to name only a few examples.Determining the runoff pathways for a study area is not evident with regard to issues of scale, resolution and water storage.Since the advent of digital elevation models (DEM), computational routines have been proposed that simplify the process in many ways and allow for at least semiautomatic procedures to delineate the pathways and drainage basins (e.g.Tarboton, 1991;Renssen and Knoop, 2000;Döll and Lehner, 2002).This has made routing runoff over a terrain a standard procedure in geographic information system applications.
Generally, all methods based on DEMs search the runoff directions using a form of steepest gradient method following the gravitational force.The level of sophistication of the runoff models is in most cases limited by available data, computational costs and specific requirements of a given situation.When specific data are available for the region of interest, different types of refinements can be used to improve the model.Examples include the "stream burning" technique (Maidment, 1996;Renssen and Knoop, 2000), which can be used when data on observed river networks is available.Known locations of river networks are deepened in the original DEM to ensure runoff is collected in those paths.Information on human-built dams and reservoirs can be included to improve the predictions (Coe, 2000) and "ridge fencing" has recently been proposed as a method to include information on known high points and divides to separate drainage Published by Copernicus Publications on behalf of the European Geosciences Union.

H. Goelzer et al.: A dynamic continental runoff routing model
basins (Masutomi et al., 2009).Apart from the runoff direction, the temporal characteristics of runoff can be captured using network response functions (Gong et al., 2009), which include the effect of local short-term water storage.
All of the described refinements depend on data availability and to a lesser extent on the possibility to compare results to present-day observations in the form of stream gauge stations.For paleo-applications, the focus of the present work, such specific data (e.g. on past river locations) are sparse or simply not available, which simplifies the task but also introduces uncertainties in the runoff predictions.In addition, details of the temporal characteristics of the routing network are not crucial to resolve.Water travel time through the network, typically of the order of days, may safely be ignored given the long response time scales of interest.
The runoff model described in this paper was conceived and tested for one specific paleo-application: the decay of the Northern Hemisphere ice sheets during the last deglaciation and the Holocene (21-0 kyr BP).Aside from the associated large changes in topography and albedo, their freshwater discharge probably had strong effects on the ocean circulation and on the global temperature evolution during this period.Paleoclimatic evidence suggests that rapid climate changes during the glacial cycles can be associated with reorganisations of the ocean-atmosphere system (e.g.Broecker et al., 1985;Clark et al., 2002).The location of freshwater entering the ocean is thereby thought to have a strong control on the evolution of the ocean circulation (e.g.Broecker et al., 1990).To study possible changes in the ocean circulation, it is therefore crucial to track both the intensity and location of glacial melt water runoff and of iceberg calving.
While the routing of freshwater for a given topography is a standard procedure, and an example exists for the time period of the last deglaciation (Tarasov and Peltier, 2005), our model is specifically aimed for and ready to be coupled to an ocean model.We therefore address the problems that arise from the combination of a temporally changing topography (as a function of evolving ice thickness and isostatic adjustment) with an ocean model grid that uses a fixed landsea mask.This represents an important step towards coupled ice-sheet-ocean model simulations.At present, ocean models that allow for a temporally variable bathymetry and land-sea mask are not readily applied on a global scale and time scales of glacial-interglacial transitions, due to computational limitations.In our application, we therefore deal with a constant ocean model grid in combination with a variable land-sea mask of the ice sheet model.The mismatch between ice sheet and ocean model grids has to be taken into account, and the changing topography makes it necessary to regularly update the routing matrix automatically and in a computationally efficient way.

Model setup
The continental runoff model (CRM) described in this paper is a new module to an existing Northern Hemisphere ice sheet model (NHISM, Zweck and Huybrechts, 2005) and serves as interface to an ocean model (CLIO, Goosse and Fichefet, 1999).Before describing the CRM in detail, we first give a brief overview of the ocean model grid and the ice sheet model.Although developed specifically for the given model setup as interface between NHISM and CLIO, the fundamental formulation of CRM could be applied to other combinations of ice sheet and ocean models and for other applications, where the routing matrix has to be updated automatically.

The Northern Hemisphere ice sheet model
The ice sheet model used for the present study merely serves to illustrate the use of the continental runoff model, so only a short description is given here.NHISM, a threedimensional thermomechanical ice dynamic model consist of three main components which respectively describe the ice flow, the solid Earth response and the mass balance at the ice-atmosphere interface (Huybrechts and T'siobbel, 1995;Huybrechts and de Wolde, 1999;Huybrechts, 2002).In many ways the model is similar to other large-scale ice sheet models (e.g.Ritz et al., 1996;Greve, 1997) based on the shallow ice approximation, which drastically simplifies the calculation of ice velocities.NHISM only considers grounded ice and therefore does not account for ice shelf dynamics.Nevertheless, marine ice dynamics is included with a parameterization for marine calving that allows determining the extent of ice grounded below sea level (Zweck and Huybrechts, 2003).Freshwater fluxes from the ice sheet consist of spatially resolved surface ablation, calving at the margin and basal melting.The isostatic adjustment module included in NHISM accounts for changes of the Earth surface due to the evolving weight of the ice.It principally consists of two separate layers of the Earth's mantle, an elastic plate (lithosphere) and a viscously deforming layer with a characteristic adjustment time scale τ underneath it (astenosphere).The lithospheric treatment accounts for a flexural rigidity D, meaning that aside from local isostasy, contributions from remote locations are taken into account as well (Le Meur and Huybrechts, 2001).Note that limitations of this type of isostatic model have been reported when applied for glacial cycle modelling of Eurasia (van den Berg et al., 2008).

The ocean model grid
In the following we briefly describe the horizontal grid of the ocean model CLIO (Goosse and Fichefet, 1999), which is the only aspect of relevance for the coupling with the runoff model.CLIO is the three-dimensional ocean-sea ice module of the Earth system model of intermediate complexity   LOVECLIM (Goosse et al., 2010).In order to avoid singularities at the poles of the coordinate system, two spherical subgrids are used in this model, for which the poles are not located in the ocean domain.The first one is based on classical longitude-latitude coordinates.It covers the Southern Ocean, the Pacific Ocean, the Indian Ocean and the South Atlantic.The remaining parts of the ocean (i.e. the North Atlantic and the Arctic) are represented on the second spherical subgrid, which has its poles located at the equator, the "North Pole" in the Pacific and the "South Pole" in the Indian Ocean.The two subgrids of both 3 ×3 • horizontal resolution are connected in the equatorial Atlantic where there is a correspondence between the meridians of the South Atlantic on one grid and the parallel in the North Atlantic on the other grid.In the model, the two subgrids are considered as a single coordinate system using metric coefficients, which are given by Deleersnijder et al. (1993).

Continental runoff model
For a given surface topography, the continental runoff model (CRM) calculates the runoff pathways and connects each continental point with an oceanic grid point to which freshwater fluxes can be delivered.The basic assumptions are that water does not evaporate nor infiltrate in the underlying layers and that enough water is available to fill up continental depressions.We furthermore ignore temporal aspects of the routing process assuming an instantaneous redirection of freshwater fluxes, which is a good approximation for the time scales under consideration here.Topographic changes are in this study calculated by NHISM and consist of both changes in ice thickness and isostatic adjustment of the Earth's crust.An overview of the model components and processing steps is given in Fig. 1 and is described in the following subsections.

Updating the relief map
To properly calculate continental drainage paths it is decisive to have a high-resolution topographic map, resolving as many features of the surface terrain as possible.Furthermore, the model domain needs to be sufficiently large to allow water to flow to the appropriate river mouth, at possibly long distances over the continents before reaching the ocean.Consequently, we decided to combine topographic changes as calculated by NHISM (50 × 50 km resolution) with a hydrologically correct present-day digital elevation model (DEM) at 1 × 1 km resolution, called HYDRO1k (http://gcmd.nasa.gov/records/GCMD HYDRO1k.html).The dataset does not contain bedrock information of, e.g. the Great Lakes and of ocean bathymetry.It was therefore combined with a 1 resolution DEM with global coverage (ETOPO1; Amante and Eakins, 2009).
To combine the HYDRO1k map with topographic changes from NHISM they both need to be transformed into a common grid.Therefore, HYDRO1k was transformed to a Polar Stereographic projection with standard parallel at 60 • N (hereafter called HYDRO1kSP) by means of 2-dimensional Lagrange polynomials.Subsequently, NHISM output was interpolated onto this new grid, which has a resolution of 6.25 × 6.25 km (Fig. 2).While the NHISM grid is limited to the region where ice was present during past glacial periods, the HYDRO1kSP grid extends further south to include potential continental routing pathways.Grid points where ice is present replace HYDRO1kSP values with the absolute surface elevation given by NHISM.For ice-free grid points within the NHISM grid, isostatic bedrock changes are added to the HYDRO1kSP topography.In order to avoid discontinuities, isostatic bedrock changes at the edge of the NHISM grid are linearly phased out to zero at the map boundary of HYDRO1kSP.

Constructing a region mask
Since CRM is constructed as interface between NHISM and CLIO, continental drainage areas have to be grouped and connected to grid boxes of that specific ocean model.While the drainage is calculated for a temporally changing landsea mask on the HYDRO1kSP grid, the runoff has to be mapped to the much coarser and fixed CLIO grid.Consequently, coastal points on the HYDRO1kSP grid may not belong to any of the oceanic CLIO grid boxes.Another problem arises from inland grid points below sea level, e.g. in the Caspian Sea, which have to be connected to the CLIO grid as well.We therefore follow the approach of Tarasov and Peltier (2005) to introduce a depth threshold and search the routing end point of any continental grid point 600 m or deeper below current sea level.Furthermore, we increase the ocean mask by extending it radially around each CLIO grid point.To connect routing end points that are still outside the extended ocean mask, e.g. in the Black Sea and Red Sea, a search algorithm is used, which locates the closest CLIO grid point.A similar complication arises from the limits of the HYDRO1kSP domain.Continental regions that are intersected by the model domain boundary also have to be connected to the closest ocean grid point using the same search routine.Given all the complications, a combined mask of the following regions is constructed on the HYDRO1kSP grid (Fig. 3): 1. Continental border region: combines all grid points above sea level at the edge of the HYDRO1kSP domain.
2. Continental region: combines all grid points above sea level except for R1.Regions with a depth between current sea level and the depth threshold (R3, black) close the gap between topographic coast (R2, grey) and extended ocean mask (R5, white).Map border points (R1, unmarked) and internal ocean regions, which are below the depth threshold but not part of the extended ocean mask (R4, red) are mapped to nearby ocean grid points.
3. Shallow water region: combines all grid points between sea level and the specified depth threshold.
4. Internal ocean region: combines grid points below the depth threshold outside the CLIO ocean mask.
5. Deep ocean regions: contain all grid points below the depth threshold, which belong to the CLIO ocean mask.
During interglacial conditions with a higher sea level stand, the distance between the topography-defined coast and the fixed ocean mask is larger compared to glacial conditions (Fig. 4).Regions below the depth threshold, which are not resolved on the extended CLIO ocean mask (R4, red) and border points (R1), are mapped to nearby deep ocean grid points (R5).

Calculation of depressions, lakes and water systems
Given the surface relief map and a sea level dependent region mask, the inventory of the continental runoff pathways can be pursued.As a pre-processing step, a ripple pattern is superimposed on the relief map (Fig. 5), in order to reduce the processing time, which arises from assembling a large number of small depressions into units of larger scale.The depth of these ripples of several centimetres can be neglected compared to the size of real topographic surface features.A substantial gain of model run time is achieved mainly when processing flat surface features, while the global drainage pattern remains the same compared to using the unprocessed relief map.
For each grid point, the surface gradient is followed to the lowest grid point in a depression that can be reached.If the lowest point is a deep ocean, internal ocean or border point, it qualifies as an end point to the total routing.All grid points that share the same lowest point are identified and form one of a large number of indexed "drainage basins".Individual drainage basins are grouped into "lake catchments", which comprise all drainage basins that are connected to the same lake when the relief is completely filled with water.At this point the outlet of each lake is determined, and the lake volume is computed for later use in the lake storage module.In turn, the lake catchments are grouped into "water systems" consisting of a cascade of lakes that have connected outlets.Water from the lowest lake in a water system can be drained to a specific ocean grid point, since the outlet of the lowest lake is by definition connected to an end point in the deep ocean, internal ocean or on the map boundary.All water systems with endpoints belonging to the same CLIO grid are taken together (Fig. 6a).Since NHISM calculates freshwater fluxes on a 50 × 50 km grid resolution, the centre of the NHISM grid box is used to define the entry point into the higher resolution HYDRO1kSP routing grid.

Lake storage module
The additional lake storage module is mostly a postprocessing routine that modifies the output to each CLIO grid cell according to reservoir changes on the runoff pathway.Calculating changes in water storage and consecutive redistribution on the high resolution HYDRO1kSP grid, although theoretically possible, is too time consuming.We therefore opted for a scheme that operates with the effective freshwater fluxes that reach the ocean by tracing changes in water storage and runoff connected to each CLIO grid box at a time.The routine takes as input the runoff and the potential volume, i.e. the integrated volume of all lakes in the connected water systems.Initially all lakes are assumed as completely filled and the integrated water storage equals the potential volume.Surface topography changes can then increase or decrease the potential volume, the latter of which leads to at times major drainage events.Runoff from the ice sheets either refills the water storage or is released when the storage has reached the potential volume.The time period over which water is released to the ocean cannot be derived from the model itself, but rather has to be chosen by the user.A guideline for this model parameter choice may be found in paleo-evidence of past large scale drainage events and furthermore has to be in line with the sensitivity of the ocean model and its capability to deal with large amounts of localized freshwater input (e.g.Manabe and Stouffer, 1997;Renssen et al., 2001).

Results
The first immediately available result of CRM is the delineation of drainage regions that are connected to the same ocean grid box (Fig. 6a), and which will be discussed in more detail in the next section.Changes of the surface topography over time due to changes in ice cover, sea level and isostatic adjustment lead to a reorganisation of the drainage regions.While these are interesting in their own rights, it is both more instructive and physically meaningful to group the different drainage regions into larger units, according to the ocean domains the runoff is directed to.Freshwater fluxes from the ice sheets during the last deglaciation are thought to have had a strong impact on the oceanic circulation (Broecker et al., 1985), while their location of impact is discussed as an important factor (e.g.Tarasov and Peltier, 2005).The division of drainage domains here is for diagnostic purposes only and somewhat arbitrary, but motivated by these research questions.The same grouping of drainage basins has been used to display the time series of runoff partition and lake drainage events from the decaying Northern Hemisphere ice sheets over the history of the last deglaciation and the Holocene (Fig. 7), which we give here for illustrative purposes only.It is beyond the scope of this paper, with focus on the CRM model description, to compare the location and timing against results of other models or geological reconstructions.
Lake drainage events (Fig. 7, bottom) that can exceed the volume flux from direct freshwater runoff arise in relation to major reorganisations of the routing matrix and when ice dammed lakes suddenly drain due to removal of the ice barrier.

Model validation for the present-day topography
For the present-day configuration the derived drainage network can be compared to major observed river networks, even if some ocean grid points in the model receive input from more than one river system.We have compared the results on a visual basis with maps from the "Watersheds of the World" study of the World Resource Institute (Revenga et al., 1988) and from "The Times Atlas of the World" (Times Books, 1998).To validate the performance of the model further, we have compared the results with the level one drainage basins of the HYDRO1k data set (Fig. 6).The delineation of drainage basins in CRM is overall in good agree- ment with the data.At this stage it is difficult to separate the influence of the ice sheet and isostatic models, which adds uncertainty to simulations and makes it difficult to evaluate the model during past periods.There is however a growing body of geomorphological data that can eventually be used to constrain the model (e.g.Toucanne et al., 2010).

Model sensitivity to grid resolution
The ultimate goal of developing the CRM is to use it inside the framework of a fully coupled global atmosphereocean-ice sheet model with other components competing for computational resources.It is therefore crucial to be able to run the model at low resolutions with a reasonable representation of the runoff routing on a scale that the ocean model is sensitive to.The loss of small-scale topographic features in coarser resolution DEMs can cause problems for the drainage calculations, e.g. for locations where narrow valleys are blocked or real barriers are removed because they are not resolved anymore.Similar issues arise for other topographyrelated problems (e.g.sills in the ocean and mountain peak altitude for atmospheric models) when working on coarse resolution, without one apparent solution.In the present study, starting from HYDRO1kSP as described above, lower resolution DEMs are derived by sub-sampling the high resolution DEM after smoothing with a two-dimensional Gaussian filter (with a standard deviation close to the target grid resolution).This results in four DEMs with horizontal resolutions of 6.25 km, 12.5 km, 25 km and 50 km.The computing time of CRM decreases rapidly with decreasing resolution of the DEM.Compared to the highest resolution DEM (6.25 km resolution), the calculations are a factor of 7, 35 and 70 faster for the 12.5 km, 25 km and 50 km resolution versions, respectively.The bulk of the increased computing time for higher DEM resolutions is spent in the search routine that assembles small-scale local depressions into units of larger scale.
Due to known choke points, the initial delineation for the present day configuration for lower resolution versions of the DEM (12.5 km, 25 km and 50 km) show some large-scale differences to the delineations in the HYDRO1k data set.In Europe the largest problem is that most of the Danube catchment area would drain into the North Sea instead of into the Black Sea.This is due to a very narrow part of the Danube valley at the border between Serbia and Romania, called the "Iron Gate", which is not resolved on the low resolution DEMs.Without modifications, the size of the Mississippi drainage basin would equally be too small in our coarse resolution models, due to the unresolved narrow valley of the Missouri river in North Dakota.Its runoff would be blocked and water would be redirected to the north instead of to the south.These differences to the observed delineations have been resolved by minor local editing of the DEM to achieve a more accurate present day basin divide.The changes made for the present-day configuration is promising to also hold for past configurations, given that ice thickness evolution and isostatic adjustment in our model change the terrain on a relatively larger spatial scale compared to the small compromising topographic features.However, Tarasov and Peltier (2006) report on the necessity to take into account changes in past sill elevations for an accurate deglacial drainage calculation on 50 km resolution.
For the comparison between different horizontal resolutions of the DEM, we exclude lake storage effects in the model, which guarantees that the total amount of runoff is identical in all four experiments at any time.This also makes it possible to trace the origin and destination of individual changes in the runoff time series.
The routing pattern on a local scale, i.e. on the level of individual water systems and of regions connected to individual ocean grid boxes, shows some variations between different horizontal resolutions of the DEM (not shown).In contrast, the amplitude of drainage to the large-scale ocean basins is largely similar between different resolutions (Fig. 8), except for the lowest resolution of 50 km.The small differences in freshwater flux between the other experiments (resolution 6.25 km, 12.5 km and 25 km) can be expected to have a minor influence on the large-scale ocean and climate model response in a coupled model setup.The above comparison therefore shows that a DEM resolution of 25 km is sufficient for our aim to model the freshwater routing for such an application.It would be difficult to justify a manifold increase in computing time using a higher resolution version of the DEM.

Model sensitivity to changes of the isostatic adjustment module
In order to study the sensitivity of CRM to changes in isostatic adjustment, we use the same parameter modifications in NHISM as described in Zweck and Huybrechts (2005).
The two adjusted parameters are the lithospheric rigidity D, which is related to the bending characteristics of the Earth's crust, and the decay time of isostatic adjustment τ , which is related to the viscosity of the mantle (Table 1).The investigated range of D of between 10 24 Nm (LR0.1) and 10 25 Nm (reference model, REF) represents the range of geophysically realistic values for the effective thickness of the lithosphere.This is complemented by an extreme case with D = 0 Nm (LR0), which is equivalent to local glacio-isostatic compensation.The sensitivity of CRM to changes in τ is examined over its full possible range of values, from instantaneous isostatic adjustment (τ = 0 ka, LT0k) to adjustment time scales of 1000, 3000 and 10 000 yr (LT1k, REF, LT10k) and no adjustment at all (τ infinite, LTinf).Similar to the procedure for resolution experiments we ignore lake storage contributions for this comparison of different isostatic model parameters, and, building on the findings above, we use a lower model resolution (25 km) in order to save computing time.Since the ice sheet evolution now differs between experiments, which complicates the comparison we prescribe a spatially constant freshwater input over the whole model domain.This is crucial for the interpretation of changes in the routing direction independent from the effect of varying local freshwater availability (Fig. 9).
The case of no adjustment (LTinf, Fig. 9g) shows changes in the routing direction as a function of ice thickness changes only.Note however that the resulting Northern Hemisphere ice sheets in this case are higher and larger compared to the reference case.The suppressed downward response of the bedrock makes the ice sheets grow higher by strengthening the elevation-temperature feedback.Transitions to the final interglacial configuration consequently appear later in LTinf compared to the case of instantaneous adjustment (LT0k, Fig. 9d) and tentatively, subsequently earlier in LT10k (Fig. 9f), the reference case (REF, Fig. 9a) and LT1k (Fig. 9e).Changes in lithospheric rigidity D (compare Fig. 9a, b and c) have less importance for the routing than changes in the astenosphere adjustment time scale τ .The overall picture that emerges is that changes in the isostatic adjustment module have a large influence on the resulting routing time series.This could have been expected due to the large range of parameter modifications especially for the adjustment time scale τ .The present-day routing partition between the oceanic basins is however very similar between the different cases, because the DEM is relaxed towards the present-day configuration.
In order to increase the robustness of the test on resolution dependence presented above, we have covered a wider range of possible deglaciation chronologies by repeating the test for two modified isostatic configurations (LR0.1 and LT10k).The results confirm independence of the large scale routing pattern and magnitude as a function of resolution reduction to at least 25 km.In fact, we find similar large-scale drainage

Conclusion and discussion
A continental runoff routing model has been presented, which dynamically updates the freshwater routing matrix in response to topographic changes and takes lake storage effects into account in offline mode.The model output was compared against present-day observed river networks and applied for illustrative purposes in the framework of simulating the last deglaciation of the Northern Hemisphere ice sheets.
It was shown that the large-scale ocean basins and timing of major freshwater forcing periods is relatively stable between versions of different DEM resolution.For the range of DEM resolutions under scrutiny, the lowest resolution sufficient for calculating the large-scale runoff pathways necessary to determine the routing of freshwater fluxes to an ocean model is 25 km.Changes in the isostatic adjustment module of the ice sheet model have a strong impact on the timing and direction of runoff routing.This is especially true for modifications of the lithosphere adjustment time scale τ .
Results from the CRM are subject to the following confinements and limitations.The routing matrix is calculated under the assumption that continental depressions are filled at any time.Furthermore, ground water infiltration, and evaporation are not accounted for.Surface erosion and sediment dams are ignored and the drainage model does not allow for sub-glacial drainage.As opposed to other routing algorithms that include dynamic lake storage (Tarasov andPeltier, 2005, 2006), our lake module operates in post-processing mode only.This introduces a source of error to the computed elevation changes since the water load is missing in the isostatic calculation.At the same time lake calving cannot be included.Finally, the low resolution versions of the DEM (12.5 km, 25 km and 50 km) do not account for small-scale topographic features necessary to resolve some known choke points, requiring additional editing (see validation section).Modifications of the DEM to match the observed presentday river networks are however promising to hold also for past configurations, since ice thickness changes and isostatic adjustment modify the DEM over relatively large distances.Note that in the present standalone analysis only freshwater fluxes from the ice sheets are routed by CRM, but precipitation over land is not accounted for.The derived routing  matrix could in principle be used in a fully coupled system to also route the runoff from precipitation over ice-free land.This would add another source of spatial and temporal variability to the freshwater forcing of the ocean and has the potential to increase the importance of routing direction changes.
As a further development for coupled ice sheet-CRM experiments, changes in water loading could be submitted to the ice sheet model to be included as an additional weight contribution in the isostatic adjustment module.When ocean models with variable bathymetry become available for largescale and long-term applications, the presented approach can be promptly adjusted to deal with dynamic land-sea mask changes. 19

Figure 1 :
Figure 1: Overview of the Continental Runoff Model and of exchanges with the ice sheet model (NHISM) and the ocean model (CLIO).The processing chain on the right hand side displays how freshwater fluxes from the ice sheets are routed to appropriate oceanic grid cells.

Fig. 1 .
Fig. 1.Overview of the Continental Runoff Model and of exchanges with the ice sheet model (NHISM) and the ocean model (CLIO).The processing chain on the right hand side displays how freshwater fluxes from the ice sheets are routed to appropriate oceanic grid cells.

Fig. 3 .
Fig. 3. Schematic delineation of different topographic regions R1 to R5. Depressions are marked d1 to d7 and lakes L1 to L5.All grid points connected to the same lowest point of a depression (example d2) form a drainage basin (red).Drainage basins connected to the same lake (example L2) form a lake catchment (green).See text for further details.

Figure 4 :Fig. 4 .
Figure 4: Ocean mask and topographic regions for the LGM (a) and the present day (b).2Regions with a depth between current sea level and the depth threshold (R3, black) close the 3 gap between topographic coast (R2, grey) and extended ocean mask (R5, white).Map border 4 points (R1, unmarked) and internal ocean regions, which are below the depth threshold but 5 not part of the extended ocean mask (R4, red) are mapped to nearby ocean grid points.6 7

Figure 5 :
Figure 5: Ripple pattern overlaid on the raw topography relief map to decrease computational costs in the solver in flat areas.The ripple depth of several centimetres is negligible compared to real topographic surface features.

Fig. 5 .
Fig. 5. Ripple pattern overlaid on the raw topography relief map to decrease computational costs in the solver in flat areas.The ripple depth of several centimetres is negligible compared to real topographic surface features.

Fig. 6 .
Fig. 6.Drainage basins for the present-day topography derived with CRM (a), where each colour represents a region that drains into the same ocean grid box (shown in grey over the ocean).(b) Level one drainage delineation from the HYDRO1k data set.

Figure 7 :Fig. 7 .
Figure 7: Freshwater contributions to different ocean basins (top) and lake drainage events 2 (bottom) for a transient deglaciation experiment.3 4 Fig. 7. Freshwater contributions to different ocean basins (top) and lake drainage events (bottom) for a transient deglaciation experiment.

Figure 9 :
Figure 9: Sensitivity of runoff partition between major oceanic basins to changes in the 2 isostatic module of the ice sheet model.a) Reference run [REF], changes of the lithospheric 3 rigidity LR0.1 (b), LR0 (c) and changes of adjustment time scale LT0k (d), LT1k (e), LT10k 4 (f), LTinf (g).Parameter values are given in Table 1.The freshwater input was set spatially 5 constant over the entire model domain in order to remove the effect of local freshwater 6 availability.See text for details. 7

Fig. 9 .
Fig. 9. Sensitivity of runoff partition between major oceanic basins to changes in the isostatic module of the ice sheet model.(a) Reference run [REF], changes of the lithospheric rigidity LR0.1 (b), LR0 (c) and changes of adjustment time scale LT0k (d), LT1k (e), LT10k (f), LTinf (g).Parameter values are given inTable 1.The freshwater input was set spatially constant over the entire model domain in order to remove the effect of local freshwater availability.See text for details.