Generating synthetic fjord bathymetry for coastal Greenland

. Bed topography is a critical boundary for the numerical modelling of ice sheets and ice-ocean interactions. A persistent issue with existing topography products for the bed of the Greenland Ice Sheet and surrounding sea ﬂoor is the poor representation of coastal bathymetry, especially in regions of ﬂoating ice and near the grounding line. Sparse data coverage, and the resultant coarse resolution at the ice/ocean boundary, poses issues in our ability to model ice ﬂow advance and retreat from the present position. In addition, as fjord bathymetry is known to exert strong control on ocean circulation and ice-ocean 5 forcing, the lack of bed data leads to an inability to model these processes adequately. Since the release of the last complete Greenland bed topography-bathymetry product, new observational bathymetry data have become available. These data can be used to constrain bathymetry, but many fjords remain completely unsampled and therefore poorly resolved. Here, as part of the development of the next generation of Greenland bed topography products, we present a new method for constraining the bathymetry of fjord systems in regions where data coverage is sparse. For these cases, we generate synthetic fjord geometries 10 using a method conditioned by surveys of terrestrial glacial valleys as well as existing sinuous feature interpolation schemes. Our approach enables the capture of the general bathymetry proﬁle of a fjord in North West Greenland close to Cape York, when compared to observational data. We validate our synthetic approach by demonstrating reduced over-estimation of depths compared to past attempts to constrain fjord bathymetry. We also present an analysis of the spectral characteristics of fjord centrelines using recently acquired bathymetric observations, demonstrating how a stochastic model of fjord bathymetry could 15 be parameterised and used to create different realisations.

The predefined distance interval is a simple finite difference parameter, and is picked so that there are enough nodes on a path to resolve it -this would be the same however we implemented the above. The ideal (but unattainable) value would be zero.
The relationship between child and parent edges is chosen to be small enough allow the path to turn quickly enough to follow the channel, and large enough so that the minimum radius of curvature is of the order of one channel width. The ideal value would be 2*pi, and restricting it simply anticipates the expectation that paths containing many loops will ultimately be rejected. In other words, this keeps paths close to the locally optimal direction, but allows some latitude so that branches can form.
The angle between new paths (pi/24) is another finite difference parameter -we cannot consider the continuum of angles spanning (-pi/6,pi/6). Again the ideal value would be zero and we reduce it progressively until the same (or similar enough) set of paths are generated.
The values |x-xi|, |y-yi| < 16 km are chosen to identify similar paths. The first two are chosen such that one of a set of seeds is identified within a small number of generations, if the paths that start from from them appear to be converging.
The condition |a-ak| < pi/8 is also chosen to identify similar paths. The angle must be greater than zero to allow branches to form from paths which are identical up to the point of the branch, and persist provided they arrive at distinct end point. A smaller angle reduces the number of branches in play at any one generation.
We have added additional text to section 3.1 to add further clarity. 5) Pages 7,8: I think that the demonstration of centerline picking would greatly benefit from using a path that was not a perfect straight line. This part is a bit hard to follow, and I think a curved path would give the reader a more intuitive sense about how it really works. I would suggest in particular that the geometry used for Figure 6 be the same as used for Figure 5. I also wonder: in looking at Figure 4b, wouldn't it be a lot simpler to just follow the path of maximum distance?
In line with the suggestions of both reviewers 1 and 2, to make the methodology as clear as possible, we have amended figure 5 so that it uses the same geometry as presented in figure  6. The algorithm is based on following the path of maximum distance, and all the additional complexity comes about because there are many possible start and end points. A simple approach might be to compute every optimal path between every possible start and end, and then choose among them somehow. For example, imagine two possible start points (seeds) and one possible end point. There are two optimal paths, and we need some way to decide whether to keep both, or just one (e.g if they are at the start of the same fjord). Our method is really just attempting to do this choosing 'on the fly', by discarding paths that start at nearby points to a 'better' path if they appear similar otherwise (have the same centroid and direction). By 'better' we mean that a path is closer to the centre of the channel (as measured by the maximum distance) on average along its length.
We have added additional text to section 3.1 to add further clarity.
6) Page 14, line 26: It is inaccurate to say that this is the first time such a methodology has been applied to fjords since subglacial channels are geometrically identical to fjords.
We agree with the reviewer and have amended the top of the discussion. The discussion now begins "Channel elevation point meshes have been implemented in different research fields including hydrology (Merwade et al., 2005(Merwade et al., , 2008, and glaciology (Goff et al., 2014). This study provides a key addition, which addresses sparse data availability with…" Response to review #2: Anonymous 1. I agree with John Goff that the description of the centreline method could be made clearer if Figures 5 & 6 were comparable, and if Fig.5 included a bend. Could you have Fig. 5a as a straight fjord to clearly demonstrate the algorithm, and then a 5b with a simple bend? Then you could probably leave Fig. 6 as it is.
We agree with both reviewers 1 and 2 on this point and have amended Fig. 5 so that (5a) now represents that which was Fig. 5(b) in the original manuscript. The new 5(b) now uses the same geometry as presented in Fig. 6.
2. I found it hard to compare Figs. 12 & 14 because of the different sizes of the plots and axis scales. I know they aren't displaying exactly the same thing, but it would help the visual analysis if the x and y axes were equivalent in scale.
We agree with the reviewer. This is a good point and we have now changed the width of Figure 14 (the high frequency stochastic model) to better match Fig. 12 (centreline elevation profiles). When doing this we realised that our stochastic model amplitudes were not correctly normalised, so we have re-done the plot 14 with correctly normalised amplitudes. The new, correctly normalised amplitudes appear significantly smaller than the real fjord elevation range (since they just deal with the higher frequency perturbation). In order to make this point clear, and to better relate Figs 12 and 14 we have modified caption 14 as follows: Figure 14: Two different realisations for the stochastic model for high frequency perturbations to the synthetic fjord elevation profiles. The model uses the parametric fit in Fig. 13 to generate the profiles, and is statistically consistent with the OBS1516 bathymetric profiles (green lines in Fig. 12). The overall trend of the fjord bathymetry and lower frequency oscillations (corresponding to wavelengths ~14 km or greater) are not synthetically generated and explains why the amplitude of the modelled elevation) is significantly less than the bathymetric observations in Fig 12. 9. page 7, l. 30-34: It is clear that some of the potential paths have to be removed along the track, but I'll have to admit I fail to understand the explanation why step 4 does this job the way it is supposed to be. Any chance to put this into some simple descriptive words illustrating the idea and the reasoning?
The lack of transparency of this step was also highlighted by reviewer 1, to whom we responded with the following: The referee is correct to note that these parameters are somewhat arbitrary. Ultimately, they are chosen so that at least one path is found for each fjord, but not many more. We plan to formalise this process for the final DEM product Essentially, what we are doing here is solving a optimization problem many times to find the optimal paths between a large set of start-and end-points, then choosing a subset of the them such that there is only one path between each pair of close start and end points. We have complicated the process somewhat by attempting to carry out the path generation and selection at the same time, in order to avoid the large number of paths that we would otherwise need to consider.
The predefined distance interval is a simple finite difference parameter, and is picked so that there are enough nodes on a path to resolve it -this would be the same however we implemented the above. The ideal (but unattainable) value would be zero.
The relationship between child and parent edges is chosen to be small enough allow the path to turn quickly enough to follow the channel, and large enough so that the minimum radius of curvature is of the order of one channel width. The ideal value would be 2*pi, and restricting it simply anticipates the expectation that paths containing many loops will ultimately be rejected. In other words, this keeps paths close to the locally optimal direction, but allows some latitude so that branches can form.
The angle between new paths (pi/24) is another finite difference parameter -we cannot consider the continuum of angles spanning (-pi/6,pi/6). Again the ideal value would be zero and we reduce it progressively until the same (or similar enough) set of paths are generated.
The values |x-xi|, |y-yi| < 16 km are chosen to identify similar paths. The first two are chosen such that one of a set of seeds is identified within a small number of generations, if the paths that start from from them appear to be converging.
The condition |a-ak| < pi/8 is also chosen to identify similar paths. The angle must be greater than zero to allow branches to form from paths which are identical up to the point of the branch, and persist provided they arrive at distinct end point. A smaller angle reduces the number of branches in play at any one generation.
We have added additional text to section 3.1 to add further clarity.
10. page 8, l. 7: I think the explanation for step 7 is easier to understand if the words "When considering the length of all complete paths," are removed and the sentence simply starts with "Where …".

Removed.
12. page 13, l. 22: I suggest to remove the paragraph break here. Maybe the sentence in l. [20][21][22] can be rewritten to clarify that the improvement is an improvement compared to Bed2013, with OBS1516 as a reference what the targeted truth is.
The paragraph break has been removed and this has now been rewritten: "Considering centreline profiles for all fjords, we illustrate the improvements made to the general elevation profile of each fjord ( Fig. 12) relative to those present in Bed2013, by considering the general agreement between the synthetic geometry and OBS1516." 13. page 14., l.1: I think "provided" or "presented" is better here than "illustrated" We now use presented.
We now ensure all uses of "high frequency" are hyphen-free.
17. caption to Figure 5: Even with the electronic version in front of me, the "Please refer to the online version of this article to make use of references to colour." bit fails to make sense to me.
This has been removed.
18. caption to Figure 8: I suggest to make it "along-transect" at both locations Done.
19. Figure 9: The figure shows that SynthIBCAO and SynthOBS give very similar results, which is good and demonstrates the ability of the method. However, I miss a possibility to directly compare to the OBS1516 data. Would it be possible to replace the grey area in panel a, which just indicates data coverage, by the OBS1516 elevation data? Keeping the extent / data coverage information of course (white areas remain white).
Yes -this is now implemented. Figure 11: I am sure the last sentence can be formulated in an easier way (still making the point to be made)

caption to
We have reworded this to: "Positive differences (red) occur where the subtrahend is deeper than the minuend, with negative differences (blue) occurring where the subtrahend is shallower than the minuend." 21. In general, I miss a statement (e.g. in the Summary) on how many manual steps are required during the procedure. Manual steps are mentioned at several locations in the manuscript, but I wonder how much work it would be to do this for, let's say, entire Greenland.
We have added this to the end of the discussion: "This paper provides a proof of concept routine for constructing geomorphologically realistic fjord geometry in the absence of observations. Actual implementation of the presented routine for large regions (e.g. the Greenlandic coast) would require manual intervention in so far as (i) identifying a seed elevation at the head of the channel and (ii) defining an end zone (e.g. the fjord mouth).
Step (i) could be achieved by using a nearest neighbour approach to acquire the nearest elevation to a given seed location. A solution to step (ii) could be by using an observation density grid where the end zone is identified as being a location with an observation density greater than a chosen value. In addition to this, the values necessary to prevent the development of closed circuit artefacts would have to be adapted to the width of the fjords for which the method is implemented."

Additional corrections
In addition to the changes made in response to the above comments, we have made some additional amendments.
Bathymetry data coverage: In line with the release of certain data from a project collaborator, some of the originally presented bathymetry data has had to be spatially restricted. The imposed restriction does not affect any of our results other than the reported statistics with regard to differences between Bed2013 and OBS1516 (section 4.1) and required changes to the data coverage presented in Fig. 7   is now slightly more restricted than was presented in the original manuscript. We have also corrected the description of the under-and over-estimates indicated by the colour scale bar. forcing, the lack of bed data leads to an inability to model these processes adequately. Since the release of the last complete Greenland bed topography-bathymetry product, new observational bathymetry data have become available. These data can be used to constrain bathymetry, but many fjords remain completely unsampled and therefore poorly resolved. Here, as part of the development of the next generation of Greenland bed topography products, we present a new method for constraining the bathymetry of fjord systems in regions where data coverage is sparse. For these cases, we generate synthetic fjord geometries 10 using a method conditioned by surveys of terrestrial glacial valleys as well as existing sinuous feature interpolation schemes.
Our approach enables the capture of the general bathymetry profile of a fjord in North West Greenland close to Cape York, when compared to observational data. We validate our synthetic approach by demonstrating reduced over-estimation of depths compared to past attempts to constrain fjord bathymetry. We also present an analysis of the spectral characteristics of fjord centrelines using recently acquired bathymetric observations, demonstrating how a stochastic model of fjord bathymetry could 15 be parameterised and used to create different realisations.

Introduction
Bed topography provides an essential boundary for modelling ice sheet dynamics, ice-ocean interactions and fjord circulation in Greenland (e.g. Vieli and Nick, 2011;Straneo et al., 2011). This widespread need for topographic information has motivated the development of digital elevation models (DEMs) for the bed topography, which combine remote sensing measurements 20 of the subglacial bed with the surrounding land and sea floor (Bamber and Layberry, 2001;Bamber et al., 2013;Morlighem et al., 2014). Each version of the Greenland 'bedmap' has provided improvements in resolution and reliability, with the most recent product to combine bed elevations and bathymetry data being Bamber et al. (2013) (from herein referred to as Bed2013). 1 The most recent Greenland-wide topography product (Morlighem et al., 2014) provides a significant improvement over previous versions towards the ice sheet margins. The development of the RTopo-2 provides another response to the limitations of Bed2013 within fjord regions, with improvements being made by including new observational data (Schaffer et al., 2016). Despite these advances, and a substantial recent increase in the amount of observational data available (e.g. Jakobsson et al., 2012;Dowdeswell et al., 2014;Boghosian et al., 2015;Rignot et al., 2016), data coverage remains poor for many coastal regions.

5
As a consequence, fjord bathymetry has not, in general, been well represented, and non-physical discontinuities between land and ocean edges are apparent. In particular, in Bed2013 physically unrealistic morphologies arise at lateral boundaries of fjord mouths, as demonstrated by examples from the Greenland coast in Fig. 1.
To address these issues, the international research community has responded by collecting and compiling a wealth of new 10 bathymetric data (e.g. Arndt et al., 2015;Boghosian et al., 2015;Rignot et al., 2016) with many other future campaigns planned (e.g. the NASA Oceans Melting Greenland (OMG) mission). It will, however, take time for extended coverage to be achieved, and some fjord regions will likely never be surveyed due to both environmental and logistical limitations associated with operating in ice-infested waters. There is, nonetheless, an urgent need to better understand and model the processes that affect the dynamics of marine terminating glaciers in Greenland and elsewhere, thus requiring fjord bathymetry to be better constrained 15 in DEMs.
Here, we present a new methodological framework for generating physically basedgeomorphologically realistic fjord bathymetry in regions of sparse observational data availability. To provide context for the introduction of our method, we first present a review of existing geostatistical approaches to interpolating channel features in DEMs (including hydrological systems, palaeo-20 glacial troughs and subglacial channels). In particular, we describe why these methods are ill-suited to regions where sparse observational data are available, which enables us to then demonstrate how our method provides a pragmatic solution to constraining the bathymetry of fjord systems. Our intent is that the presented approach will eventually be up-scaled to all unmapped fjords along the Greenland coast. This will significantly improve existing DEMs of bed geometry beneath and at the margins of the Greenland Ice Sheet as well as its surrounding surface topography and bathymetry. A novel feature of the method, which 25 is inspired by analogue studies of glacial troughs (Coles, 2014), is the incorporation of predefined cross-sectional channel geometry to provide a geometric structure that is physically realistic in the absence of observations, in turn providing realistic topography for applications including ice sheet modelling.

Past approaches for interpolation and integration of channel geometry in DEMs
For the purpose of integration in DEMs, fjords (Syvitski et al., 1987), river channels and glacial troughs (Batchelor and 30 Dowdeswell, 2014) can be considered as pseudo-linear channel systems that have directional flow. In the absence of adequate direct observations, the integration of anisotropic morphology is highly desirable when interpolating channel systems in DEMs. Where observations are available, there exist methods which can interpolate additional elevations of channel features (e.g. Herzfeld et al., 2011;Goff et al., 2014). However, where there are no data available, other than the known existence of a feature (discernible from remote sensing imagery), complications arise in how to accomodate the features in DEMs. In the case of Greenland, the last data product to provide a continuous bed-to-bathymetry DEM (Bed2013) used different approaches to interpolate different topographic regions. Kriging interpolation was used for the interior of Bed2013. The bathymetry was taken from the International Chart of the Arctic Ocean (v3) (Jakobsson et al., 2012) referred to as IBCAO from this point for-5 wards. The IBCAO DEM was developed from bathymetric observations using spline interpolation following Jakobsson et al. (2012). For Bed2013, triangulation (linear interpolation) was used to predict bathymetry within the fjords between the IBCAO and interior Greenland bed DEM datsets (Bamber et al., 2013), as these regions were unconstrained by observations. Using traditional isotropic interpolation approaches, such a lack of data often results in the generation of interpolated surfaces that fail to represent true channel geometry, and often appear artificially smooth. In the case of Bed2013, this problem resulted in To capture the appropriate geometry of channels, several different approaches have been developed involving geometric (e.g. 15 Goff and Nordfjord, 2004;Merwade et al., 2005), mathematical (e.g. Herzfeld et al., 2011) and mass conservation (Morlighem et al., 2014) solutions. To place our study in the context of other interpolation methods we review previous approaches with a particular focus on resolving curvilinear features (channels). Additionally, stochastic perturbations to Greenland Bed DEMs can be employed in a variety of different ice-sheet modelling contexts (cf. Durand et al., 2011;Seroussi et al., 2011;Goff et al., 2014;Sun et al., 2014). It is possible that there will be a future need for similar stochastic modeling of fjord bathymetry, and 20 we also discuss this here.

Kriging
The key issue with interpolating features for which orientation is important (e.g. channels) is the ability to incorporate direction into the method used to develop them from observations. Kriging -a method of interpolation for which the interpolated values are modeled by a Gaussian process -is often employed to create continuous surfaces from point data (e.g. Hock and Jensen,25 1999; Bamber and Layberry, 2001;Le Brocq et al., 2010;Bamber et al., 2013). The approach accounts for the statistical properties of observations within a local search neighbourhood using a variogram function (Deutsch and Journel, 1998) and, u. Using this it is possible to incorporate various types of anisotropy within the basic framework (Merwade et al., 2006).
However, the method only holds when applied over regions sharing the same overall statistical properties whether that be, for example, the same geologic rock type or the same directional bias. When anisotropy is defined relative to a fixed Cartesian 30 coordinate system, and where data are sparse, kriging is impractical for sinuous features with constantly varying direction such as channels (see also Fadlelmula F. et al., 2016). Specifically, dividing a region into areas of shared anisotropy (thus satisfying the assumption of stationarity within a search window) that are data sparse prevents the adequate population of the variogram with which to statistically model the region.

2.2 Channel coordinate transformations
To enable interpolation across channel widths, one approach uses cross-sectional profiles, but to do so, typical channel sinuousities present a problem. As an intermediary step to interpolating sinuous channels in DEMs, several approaches have been developed to transform the coordinate system of a given channel -moving from Cartesian coordinate space to channel coordinate space -enabling removal of complex sinuosity and the creation of an artificially straight channel (cf. Goff and 5 Nordfjord, 2004;Merwade et al., 2005). Channel space (sometimes denoted as s, n in the literature) differs from Cartesian space in that locations are transformed relative to their distance along the channel (s) and perpendicular to the centreline (n).
Observations within channel space -a now straightened channel -can be locally interpolated by considering a single direction as opposed to a continuously changing one. The interpolated channel can then be transformed back to Cartesian space.
The approach breaks down, however, where multiple channels merge together at confluences. Furthermore, in the absence of 10 sufficient observations, such an approach cannot be used alone to predict along-channel geometry without additional interpolation. For example, manual digitisation has been applied to individual channels to assist in the development of a realistic bed topography for Thwaites Glacier, West Antarctica (Goff et al., 2014). Additionally, channel straightening through coordinate transformation becomes difficult where channels manifest high levels of sinuosity or sharp changes in direction (Goff and Nordfjord, 2004).

Mathematical morphology
Further issues with regard to maintaining morphological characteristics of channels, in particular ensuring known depths are honoured, are apparent in large and low resolution datasets particularly where interpolation methods are applied (Herzfeld et al., 2011). Where resultant data products are to be used in modelling studies, honouring known maximum depths is key as incorrect values can adversely affect results -especially with regard to maximum and minimum elevations (Herzfeld et al., 2011). To directions provides a succinct approach to avoid the constraints of regular gridding, which mask channel structures especially at lower resolutions. However, observations are required to identify channels and application of the 'mathematical morphology' approach becomes complicated in the case of multiple interconnected dendritic type networks.

Mass conservation
Subglacial channels, which occur beneath grounded ice, are significantly easier to interpolate into DEMs than fjords as a 30 physically based mass conservation optimisation scheme can be applied (Morlighem et al., 2011(Morlighem et al., , 2014. This approach is independent of traditional geostatistical interpolation methods. Bed elevation values are calculated from ice thickness values, 4 which are derived from combination of radar sounding measurements and surface velocity observations and of course using the assumption that mass is conserved along flow. Despite such an approach being useful for subglacial channels covered by grounded ice, this approach cannot be applied for regions of open ocean or non-grounded ice as is the case for fjords and cross-shelf troughs on formerly glaciated continental shelves.

Remaining issues 5
Despite the approaches that have been developed to interpolate channels in DEMs, there are a number of recurring problems in applying these methods in the next generation of the Greenland DEM. In particular, all of the methods assume that there are at least some data from which to extend and predict the structure of a given feature. Furthermore, no method is explicitly designed to include or represent the known physical characteristics, in particular the cross-sectional geometry, of the particular type of channel system (e.g. U-shape of glacial; v-shape of fluvial), with morphological information only being extended from 10 available observations. Thus, there remains a disconnect between the presented frameworks and cases where features (1) are known to exist; (2) are assumed to conform to a structure related to the processes by which they were created (e.g. an assumed u-shape in the case of fjords where no other data are available); and (3) have no observations available to provide geometric constraints. A framework for fjord channel systems which addresses these issues, and can be applied to a large area such as the Greenland coast, must be able to:

15
-Impose morphological geometry to features of known process-origin; -Account for elevation trends along and across the channel; -Account for confluences in dendritic channel systems; -Enable repeatable application across numerous channels within dendritic systems; -Be able to deal with minimal data input (other than absolute limits e.g. minimum and maximum depths as well as spatial 20 extent).

Stochastic models
Stochastic models of bathymetry have long been employed to abyssal hill features in the deep ocean (e.g. Goff and Jordan, 1988). In such places, stochastic models are appropriate for use because the frequency power spectra of deep ocean bathymetry follows well defined parametric relationships (Bell, 1975). Specifically, the high-frequency tail of the power spectra is char- 25 acterised by power-law relationships (i.e. the Brownian regime, which can be stochastically modelled), with lower-frequency behaviour characterised by a flat spectraflat region of the power spectra (i.e. the white regime, which cannot) (Goff and Jordan, 1988). This spectral behaviour is common across other types of natural terrain and, subsequently, spectral analysis of natural terrain often focuses upon establishing the transition between high-and low-frequency behaviour, and the characterisation of the high-frequency power-law relationships (Shepard et al., 2001). To the best of our knowledge, there are no studies that have 30 5 considered the spectral analysis of fjord bathymetry.Whilst the spectral properties of mid-ocean bathymetry (Bell, 1975;Goff and Jordan, 1988) and subglacial channels (Goff et al., 2014) have been assessed, to the best of our knowledge this has not been done for fjord bathymetry. As part of this study, we use data that are available from surveyed fjords to constrain the stochastic models of the bathymetry of many Greenland fjords.

5
A flow diagram for the separate components of our method for generating physically basedgeomorphologically realistic fjord bathymetry is presented in Fig. 2. Each component is described in a separate sub-section. In Sect. 3.1 we discuss the approach taken to map the centreline of each fjord within the fjord system introduced below. Using the mapped centreline, we explain OBS1516. Consequently, we differentiate these simulations by naming them SynthIBCAO and SynthOBS respectively.
The sequential approach defined in Fig. 2 was applied to a fjord system in North West Greenland close to Cape York (see

Centreline mapping
The ability to map a given fjord where no observations are available requires the provision of a skeleton mesh, which hinges on the presence of a centreline -an imaginary line that is equidistant from the two fjord edges. Consequently, the first step in the synthesis of a given fjord's geometry requires a centreline to be defined. Approaches exist for automatic centreline 25 identification for glacier surfaces (e.g. Kienholz et al., 2014;James and Carrivick, 2016) as a means of avoiding manual digitisation. Such applications are, however, informed by the availability of a glacier surface elevation DEM. An equivalent non-geomorphologically based method includes the definition of the medial axis (cf. Blum, 1967) or topological skeleton and is frequently used in image processing and computer graphics applications (see Bai et al., 2007). Various packages are available to calculate topological skeletons (e.g. Van Der Walt et al., 2014). However, these algorithms are based purely on an input im-30 age and are sensitive to image pixel resolution. For our intended application, this can result in the development of a centreline 6 (or skeleton) with multiple branches along a single channel feature.
The centreline mapping method that we introduce allows fjord systems with multiple branches to be accounted for. Each centreline extends from a predefined seed point (or points) at the head of the fjord, ending at a predefined end zone (e.g. the fjord mouth). The centreline itself is defined by a series of points or vertices, each with a unique identifier. Fjord confluences 5 and the implementation of network structure are described in Sect. 3.4. We define the centreline as being any path between a seed and the end zone which minimises the path integrallength whilst maximising the distance of the path from the fjord walls.
This removes the issues of multiple side branches that arise using existing skeletonisation approaches. The algorithm incorporates direction and thus an aspect of evolutionary landscape process knowledge, which ensures that the centreline captures a leading order feature from the landscape it represents. Furthermore, this approach ensures that the paths and vertices are given 10 unique identifiers enabling them to be specifically referenced, which is important when defining the channel mesh (see Sect. 3.2).
The generation of centrelines would be straightforward if we knew in advance the start and end points for each. In that case, we would simply compute a path that minimises the line integral where C is the entire centreline, and L is some function that grows towards the channel edge. However, there are a large set of potential start points for every centreline, because we do not know ahead of time where any given channel should start.
There are also a large set of potential end points, for the same reason and also because there may be more than one centreline originating at any given start point, if the channel is forked. 20 Fjords were identified as channels between areas of land and ice leading towards the open ocean, identified here using the modified GIMP land classification product Morlighem et al., 2014) (see Fig. 3(a)). At the head of each fjord, multiple seeds were manually created from which to initiate a path. The end target was defined as a broad region rather than a specific point (in this case, the edge of the land classification mask). The following algorithmic steps were then 25 undertaken: 1. Using the land classification mask, we calculate the distance transform between land/ice and oceanof all locations between land/ice and ocean within the channel (d), from which the shortest distance of any location within the ocean from regions of land or ice can be identified (see Fig. 4).
2. Based on the slope of the distance-transform calculated for regions of ocean relative to land/ice land categories using distance intervalfinite segment length ∆l, chosen to resolve the path with sufficient detail, which propagate along the fjord (see Fig. 5(b)). Up to four new nodes are generated at each step, such that the angle between the newly defined edgesegment and the parent edgesegment is less than π/6∆l/r C , the angle between any pair of new edgessegments is no less than π/24∆a and the new edgesegment does not cross the fjord boundary. r C is chosen so that the minimum radius of curvature of any portion of the path is comparable to a typical channel width. The finite angle difference, ∆a, 5 like ∆l, is chosen to be small enough to describe the channels adequately. See table 1 for the values used. If we knew that no path would branch, we could generate a single new node: this more complicated procedure is adopted because we do need to consider branches. If more than four nodes are generated in this manormanner, then those with the smallest values of L = 1/d 4 (distance from the fjord edges) are selected.
3. Where, for example, a seed generates three new points, this results in the creation of three paths. Paths then increase in 10 length as more points are created, with new paths following the creation of each new point. In the example illustrated in 4. The process in step 3 alone would lead to exponential growth in the number of paths. To avoid this, paths are culled frequently (every three generations). Each path is categorised into bins (x i ,y j ,a k ), where the centroid of the path x, y 15 satisfies |x − x i |, |y − y j | < 16 km|x − x i |, |y − y j | < l bin , and the angle defined by the last edge added to the path, a, satisfies |a − a k | < π/8|a − a k | < a bin . The path with the lowest value for the path integral of L is retained from each bin, and the remainder discarded. l bin is chosen so that the number of paths originating from seeds at the head of the same fjord are reduced to one in a few generations. a bin is chosen so that branches originating from the same start point persist long enough to have distinct centroids if they follow genuine branches in the channel. See table 1 for the values 20 used.

5.
Where a path meets a boundary that is not the predefined end zone (e.g. land), the path is culled, as illustrated for branches 2.1.1 and 2.1.2 in Fig. 5(b) within the pink box. In this example, as a consequence of the boundary intercept, there is a resultant culling of both paths 2.1.1 and 2.1.2, following the removal of the parent node 2.1.
6. Once a given path reaches the target region, its length is compared to the length of all other complete paths with only the 25 shortest being retained. In Fig. 5(b), paths 1.11.1, 1.11.2 and 1.11.3 complete, the shortest where L is minimised (1.11.2) being retained and used to define the fjord centreline. 7. A centre of mass (COM) is calculated for each path. When considering the length of all complete paths, wWhere the distance between the COM of separate paths is greater than a threshold value (manually set to ∼ half of the mean channel width in an area), both paths are kept regardless of length, otherwise the shorter, where L is minimised, of the two paths 30 is retained. The use of the COM allows for separate centrelines to be defined along more complex fjord networks than by culling according to path length alone.

3.2 Fjord mesh development
For each fjord centreline, points normal to each centreline vertex were defined up to the fjord edge taken from the GIMP land mask (Fig. 6). The angle of the normal vector along which these new points were defined was calculated from its orthogonal relationship with the vector joining the neighbours of a given centreline vertex. To avoid an irregular distribution of new points in the interpolated profile at the channel edges, the points used to define the vector from which the normal was calculated 5 were sometimes selected from more distant neighbours. This was particularly pertinent at more sinuous sections of a centreline (see Fig. 6). This smoothing of the profile is adapted from Goff and Nordfjord (2004). Vertices normal to the centreline were calculated up to the mouth of the channel, at which point the fjord centreline was manually clipped.

Mesh elevation definition and cross-sectional fjord geometry
Elevations were attributed to the point mesh by first constraining the seed and fjord-edge bed elevations through association 10 with the nearest bedrock/bathymetric observation (either from ice penetrating radar where ice covered (see Bamber et al., 2013); altimetry where exposed bedrock ; or from bathymetric observations (OBS1516)). This method ensured that at the head of the fjord, three elevations were available -the two edges of the fjords (taken as the elevations at the first land locations encountered at the fjord edges) and a centreline elevation. For future applications along the Greenland coast where seed data are sparse, modelled estimates from the mass conservation optimisation scheme (Morlighem et al., 2014) 15 could be used.
Two approaches were taken to assign bed elevation values along the centreline of a given fjord. In both approaches, bed elevations were linearly interpolated between a known bed elevation at the head of the fjord (taken from OBS1516 (Fig. 3)) and a known bed elevation at the mouth of the fjord -the mouth having been manually located and consistently used for all 20 model runs. For the first run (SynthIBCAO), the bed elevation at the mouth was taken from the nearest IBCAO observation (20 km from the mouth of the fjord system depicted in Fig. 3) and was set at -803 m. This was chosen for the first simulation as until recently IBCAO provided the most extensive bathymetric dataset for Greenland and the distance from a fjord head to the nearest observation is often ∼ 10s of kilometres. For the second run (SynthOBS), the gridded bathymetric observation from OBS1516 at the same position was used (-920 m -see Fig. 3). Should high frequency stochastic perturbations wish to be 25 added along the profile (see Sect. 3.5 and 4.5), they would be applied at this stage. Bed elevations up to the termini of most glaciers in Greenland, albeit predominantly modelled, are now available (Morlighem et al., 2014). We justify the use of the OBS1516 data for defining the elevations at the head of each fjord in the presented simulations as it enables a comparison of synthetic and observational data directly, removing the need to consider uncertainties inherent of modelled elevations.

30
In the absence of large scale studies on fjord bathymetric geometry, we base our cross-sectional fjord geometry on the prior analysis of over 8000 glacially eroded valleys now exposed by interglacial ice-sheet retreat (Coles, 2014). In their study, profiles were acquired from different glacial and geological environments including valleys from the Southern Alps (New Zealand), the Pyrenees and North and South Patagonia. For the valleys the bed elevation, V d , was fitted to a power law relationship for the form where α is the form ratio (valley depth/valley top width), w is the distance along the cross-section from the centreline (the position of which corresponds to w=0), and β is the power law exponent. Best fit parameters, of α=0.20, β = 1.38 were ob-5 tained (Coles, 2014). A value of β = 2 (i.e. a parabolic relation) follows Wheeler (1984).
Equation 2 assumes that a given fjord's cross-section is symmetrical about the centreline with the centreline as the deepest point; an assumption which usually does not hold exactly. Additionally, the fjords are often seeded with edge elevation data that are significantly higher on one side of the fjord than the other. To define the cross-sectional fjord geometry, we define a 10 parabola of the form ax 2 + bx + c (Wheeler, 1984), where a, b and c are calculated based on the known elevations and relative locations of the edge and centreline points. This enables us to relax the constraint that the centreline must be the deepest point.
Thus, the parabola used to define across-fjord geometry in this study is inspired by the analyses of Coles (2014) but not a direct application of it.

15
Following the development of a complete fjord elevation mesh, a surface was then made for each fjord, with mesh point elevations being mapped to a regular grid, thus creating a continuous surface. The resultant grid was then masked using the GIMP mask, thus removing any values outside the extent of the point mesh which arise as a result of the regular gridding process. The individual fjord grids were then combined, from which a final grid of the minimum values (or maximum depths) was created. Thus, the lowest value at any location where two fjords overlap was retained (Goff and Nordfjord, 2004). As 20 a result, deeper grid values took precedence over those that were shallower. This approach coupled with the aforementioned setting of edge elevations at confluence locations avoids the creation of ridge artefacts in the final DEM. The fjord DEM was then integrated into the wider landscape DEM (Bed2013) which includes non-fjord regions. Prior to the merge, Bed2013 was masked, removing any values in the area occupied by the synthetic fjord(s).

Stochastic modelling of fjord bathymetry 25
In this section we describe spectral analysesanalysis methods used to constrain the fjord's statistical features, and the inverse methods that can be used to generate synthetic profile. Our analysis is based upon analogous analysis of abyssal hill features in the mid-ocean (Bell, 1975;Goff and Jordan, 1988), although it is simpler in the respect that fjord bathymetry is approximated as a one dimensional problem. Using the centreline mapping approach presented in Sect. 3.1, centreline points were established, with vertices on a 150 m interval for nine fjords along the west Greenland coast, selected where mean gridded 30 observations were contiguous along fjord centrelines according to OBS1516 (Fig. 7).
The lengths of the nine fjord sections in our example are constrained by the length of the shortest fjord section (∼30 km).
Elevations were extracted for each centreline vertex, providing one dimensional centreline elevation profiles. Where a centreline contained missing data at a level no greater than 20 %, a cubic interpolation routine was applied to give a continuous 5 elevation profile (Fig. 8). Prior to performing the spectral analysis, each elevation profile was linearly detrended, which acts to emphasise the overall variation of the small scale trends (Shepard et al., 2001). Each elevation profile was then transformed using the numerical Fast Fourier Transform algorithm converting it to the frequency domain (Van Der Walt et al., 2011). Power spectra for each fjord were then obtained from the square of the complex modulus, and arithmetically averaged over the nine selected fjord profiles, to create a composite power spectra. This arithmetic averaging approach is as described in Bell (1975) 10 for mid-ocean bathymetry, and enables longer wavelength features to be statistically constrained along with the higher frequency features that are repeatedly sampled.
Of interest in this study is demonstrating, in a proof-of-concept manner, how the composite power spectra for the fjords can be used to generate different one dimensional realisations of synthetic bathymetry, that are consistent with the overall 15 statistical properties. In order to generate the different realisations of bathymetry, we use the inverse Fourier transform method outlined by Tyan et al. (2009) (the sinusoidal approximation method described in Sect. 4 of their study). Their method was introduced in the context of generating one dimensional random road profiles, which is a mathematical analogue of fjord profiles. In their formulation the Fourier amplitudes of each harmonic are determined by the power spectra of the profile, with stochasticity present via the random relative phase of each harmonic. Our only modification to their method is to use a different 20 parametric form for the power spectra, which is motivated by our observed results described in Sect. 4.5 and is consistent with the generality of their method.

Results
In this section, we present differences between Bed2013 (the last Greenland bed-bathymetry combined data product following Bamber et al. (2013)), OBS1516 (recently acquired fjord bathymetry data), and synthetic fjord bathymetry developed using the 25 methods described in Sect. 3.1 -3.4. The first synthetic application is preconditioned using the nearest IBCAO bathymetric observations (SynthIBCAO) and the latter using the OBS1516 dataset (SynthOBS). The results are compared to Bed2013 and OBS1516, in the region illustrated in Fig. 3. The extent of OBS1516 relative to the area of interest is displayed in Fig. 9(a).
Spectral analysis of selected fjord profiles shown in Fig. 8 is then presented, following application of the method described in Sect. 3.5 to available fjord bathymetry data in the region illustrated in Fig. 7. Prior to investigating the improvements of the synthetic algorithm over the existing Bed2013 DEM within the area of interest proximal to Cape York ( Fig. 9(a)), we consider areas of maximum over-and underestimation of bed-elevation that are present in Bed2013 within the region covered by OBS1516. Bed2013 is a continuous DEM extending from the bed beneath the contemporary ice sheet out to the continental shelf in the ocean, with all bathymetric information derived from IBCAO.

5
IBCAO was combined with the bed elevation component of Bed2013, with triangulation used as an interpolator to provide values where IBCAO was unconstrained by data (Bamber et al., 2013). Triangulated portions of the resultant DEM were then smoothed using a 2 km window (Bamber et al., 2013). Where there was an unrealistic offset between the two surface datasets (e.g. bathymetry was higher than the glacier bed), some areas were manually dropped to force them to adhere to a subjectively more realistic profile (i.e. a fjord would be lower than the glacier bed upstream of it). The result of differencing Bed2013 from 10 the OBS1516 data set is presented in Fig. 10(a) with the frequency distribution of the differences presented in Fig. 10

Bathymetry for SynthIBCAO
The first implementation of the synthetic fjord routine, SynthIBCAO, defines the elevation at the mouth of the fjord based on 20 the nearest IBCAO bathymetric observation, with the elevation of the point at the head of the fjord taken from the OBS1516 dataset . Points normal to each centreline vertex were then calculated as described in Sect. 3.2 and 3.3. The resultant combined surface DEM with the inclusion of synthetically created fjord bathymetry is displayed in Fig. 9(c).
The SynthIBCAO channel geometry is both deeper and more concave compared to Bed2013 ( Fig. 9(b)), particularly with 25 regard to the narrower fjord regions. Based on the contour pattern ( Fig. 9(c)), these narrower fjord regions now display a deeper and more concave cross-sectional profile than was rendered in Bed2013. For the wider confluence region centred south from  Comparing the difference between Bed2013 and SynthIBCAO ( Fig. 11(a)), the latter dataset has elevations consistently lower than the former. The mean offset between the two datasets was 274 m. The changes along the narrower portions of the 12 fjords -up to ∼ 3 km from each respective fjord head (see Fig. 3(b) for mapped channel centrelines) -are relatively small (∼ 0 -50 m). Larger offsets are apparent where fjords enter the confluence region centred south from (-705, -1340) on Fig.   11(a), with a maximum offset of 547 m. The increased concavity of SynthIBCAO is well illustrated with a mean increase in depth along the confluence zone centreline of ∼ 370 m.

5
Subtracting SynthIBCAO from OBS1516 ( Fig. 11(b)) reveals a mean offset between the two datasets of around 50 m. Relatively good agreement along the first ∼ 3 km of each fjord (see Fig. 3(b) for mapped channel centrelines) is displayed, and indeed portions of the main confluence region, with differences centred at 0 m. The main region of synthetic elevation overestimation (i.e. lower than the observations) is focused at the confluence point of the fjord 1 (refer back to Fig. 3(b) for fjord numbers), the region of overestimation focused around (-709, -1344) on Fig. 11 For reference, a comparison with OBS1516 subtracted from Bed2013 for the same area of interest is drawn (Fig. 11(e)). 15 As described in Sect. 4.1, Bed2013 consistently underestimates bed elevation. However, as with SynthIBCAO, the main areas of underestimation are focused at the same locations -namely (704, -1338) and (-704, -1355) on Fig. 11(e) for which overdeepenings are likely present.

Bathymetry for SynthOBS
The second implementation of the synthetic fjord routine, SynthOBS, defines the elevation of points at both the head and mouth 20 of the fjord based on gridded elevations from OBS1516 at the same location. The resultant combined surface DEM with the inclusion of synthetically created fjord bathymetry is displayed in Fig. 9(d). SynthOBS demonstrates deeper concave geometry across the fjords compared to Bed2013. The changing relief of the banks of the synthetic fjords are steeper than those rendered in the original DEM ( Fig. 9(b)). Between fjords, there are also changes in the elevations of the ridges such as at (-706, -1330) on Fig. 9(d). The differences between the synthetic and the original DEMs are further quantified by the difference plot illustrated 25 in Fig. 11(c).
The SynthOBS surface is generally lower than Bed2013 (Fig. 11(c)), with a mean offset between the two datasets of 316 m. The only locations where SynthOBS was higher than Bed2013 were at the edges of the fjords within ∼ 3 km from each respective fjord head (see Fig. 3(b) for mapped channel centrelines). This possibly highlights overly smoothed sections of 30 Bed2013 (where it was combined with the IBCAO dataset -see Bamber et al. (2013)). As with the SynthIBCAO approach, the largest offsets are apparent as fjords enter the confluence region centred south from (-705, -1340) on Fig. 11(c), with a mean offset in this region of ∼ 400 m.
Subtracting SynthOBS from OBS1516 (Fig. 11(d)), the mean offset between the two datasets is ∼ -3 m. The spatial pattern is very similar to that described for SynthIBCAO with the same regions of under-and over-estimation being equally apparent.

5
Overall, the profile represented by SynthOBS is closer to OBS1516, as would be expected considering the elevation profile of each fjord was calculated between fjord head and mouth observations from the OBS1516 dataset.

Centreline profile changes: Bed2013, OBS1516, SynthIBCAO and SynthOBS
Considering centreline profiles for all fjords, the improved general elevation profile of each fjord using the synthetic approaches -SynthIBCAO and SynthOBS -relative to OBS1516, as well as the underestimation of elevation in Bed2013, are illustrated 10 in Fig. 12.Considering centreline profiles for all fjords, we illustrate the improvements made to the general elevation profile of each fjord (Fig. 12) relative to those present in Bed2013, by considering the general agreement between the synthetic geometry and OBS1516. The synthetic realisations underestimate observed bathymetric elevation to a much lesser extent than Bed2013, capturing the generally sloping profile of OBS1516. The good agreement (approximately ± 50 m) of synthetic-observed values along the first ∼ 3 km of each fjord -in particular for fjords 1,2 and 5 -implies the presence of approximately linear 15 profiles. Larger differences -indicative of where the synthetic approach performs less well -occur from ∼ 4 km along each centreline, which relates to the confluence region of the individual fjords (refer to Fig. 3). Higher frequency features (along track peaks and troughs likely relating to sills and overdeepenings) are not captured using the presented synthetic fjord bathymetry generation approaches.

20
Following Sect. 3.5, we now consider the spectral characteristics of the fjord bathymetry along the centrelines of the nine fjords illustrated in Fig. 7, using the OBS1516 data. A log-log plot for the mean power spectra, S(k) where k is the wave number (linear spatial frequency), is illustratedpresented in Fig. 13 (blue crosses). The power spectra exhibits an approximate powerlaw relationship at higher frequencies (corresponding to a linear relationship in log-log space), and an approximate flattening at lower frequencies. A parametric model which captures this frequency transition is: where k 0 represents the approximate transition frequency between the high and low frequency regimes, α is the exponent for the high frequency tail (for k >> k 0 , S(k) ∝ k −α ), and F 0 acts as a normalisation constant. The parametric model (equation 3) is a generalisation of the model for the power spectra of ocean bathymetry in Bell (1975), which assumes α=2. In general, different types of natural terrain can exhibit a range of spectral exponents (Goff and Jordan, 1988;Shepard et al., 1995Shepard et al., , 2001, 30 14 and our parametric model is representative of this. The parametric best fit values were obtained using a non-linear least squares solver, and correspond to F 0 =17.6 m 2 km −1 , k 0 = 0.069 km −1 and α=1.74 (Fig. 13, red solid line). The transition frequency, k 0 = 0.069 km −1 , corresponds to a transition wavelength of 14.5 km. This compares with a transition spatial frequency k 0 =0.025 km −1 , and a transition wavelength 40 km, 5 for abyssal hill features in the mid ocean in Bell (1975). Figure 14 shows two different realisations of synthetic fjord bathymetry using the parametric fit to the power spectra in Fig. 14, and the stochastic inverse Fourier transform procedure describe in Sect. 3.5. The horizontal spacing of the synthetically generated profiles is set to be the same as the bathymetric data (0.2 km). If we draw a comparison between the stochastic model 10 of synthetic fjord centreline profiles and the OBS1516 profiles (Fig. 12), it is clear that the synthetic profiles do not contain the lowest-frequency oscillations (wavelengths ∼ 15 km or greater). This is consistent with the general flattening of the fjord power spectra at low frequencies. However, oscillations on a length scale ∼ 5 km (typical of sills and overdeepenings) are present in the synthetic profiles, albeit the specific locations of such features in these profiles are random.

Discussion
Despite the concept of channel elevation point meshes not being new in hydrology (Merwade et al., 2005(Merwade et al., , 2008, and in some glacial studies (Goff et al., 2014), this study is the first time such an approach has been applied to large fjord systems. A key addition presented in this study, which addresses sparse data availability, isChannel elevation point meshes have been implemented in different research fields including hydrology (Merwade et al., 2005(Merwade et al., , 2008, and glaciology (Goff et al., 2014). This 20 study provides a key addition, which addresses sparse data availability with the introduction of parabolic cross-sectional form along each profile that is characteristic of glacial fjords. In the absence of data, continuous DEM surfaces are developed using interpolation procedures. The specific values assigned to regions lacking observations are thus entirely dependent on the interpolation routine applied and the presented approach provides a physically basedgeomorphologically realistic estimate of elevations in these regions. The introduction of the artificial mesh removes the need to apply a traditional interpolation routine 25 over a large region, instead providing an idealised mesh to constrain regions known to be fjords. The method presented must, however, be semi-informed by data. The minimum elevations that are required are the fjord bank edges (i.e. topographic elevation at the land/ocean interface according to a land mask e.g. GIMP) -which in general can be different from one another -as well as the elevation of the assumed centreline. The deepest point along the channel is constrained by the quadratic fit. In the case of Greenland, for which this method has been developed, ice-free edge observations are widely available (e.g. 30 Howat Korsgaard et al., 2016). Equally, observations at the head of the fjord can be taken from bed elevations inferred from mass conservation (Morlighem et al., 2014), or, in some regions, radar observations (e.g. Gogineni et al., 2001).
Finally, observations for the fjord mouth could be taken from datasets including IBCAO or others (e.g. Schjøth et al., 2012; the fjord mouth itself, which using the presented approach may result in further under-or over-estimation of a given fjord centreline elevation profile. The synthetic approaches -SynthIBCAO and SynthOBS as presented in Sect. 4.2 and 4.3 respectively -represent two sit-5 uations that would be encountered when applying the method, as part of wider Greenland DEM development, to fjords around Greenland. By informing the mouth elevation on IBCAO observational data at a distance of ∼ 20 km from the mouth, the impact of using distant bathymetric observation is exemplified. Equally, as many fjords have at least some information following various recent campaigns (including Schjøth et al., 2012;Dowdeswell et al., 2014;Arndt et al., 2015;Rignot et al., 2016), the use of observational data to constrain the algorithm is illustrated by SynthOBS. The application of these two synthetic 10 approaches has provided bathymetry more representitive of the observed elevation profiles (OBS1516) of fjords within the area of interest (Fig. 12). Within this region, topographic features, such as sills and overdeepenings, captured within OBS1516 occur. It is not possible to predict oscillatory features such as these with the geometrically flat surfaces assumed by our basic algorithm.
In these examples, the overdeeping features and sills have a length scale ∼ 5 km which is less than the transition wavelength for the fjord power spectrum of ∼ 15 km (Fig. 13). The transition wavelength provides an approximate upper bound 15 upon the length scale of features which could be modelled using our stochastic framework. Subsequently if we integrated the analysis here with the stochastic model, the overall statistics of the overdeeping features would be reasonably well represented, but their geographical locations would not.
With regard to confluences, and following Goff and Nordfjord (2004), where single channel elevation surfaces overlap we 20 accept the maximum depth. This introduces a hierarchical element to surface prediction, whereby deeper channels are favoured over shallower ones. However, as the approach is based solely on topography (not rock type or age as such information are rarely available), this introduces a limitation that cannot easily be resolved in light of such sparse observations. We suggest that, in the absence of data, use of the deepest value is preferable over shallower values, due to the overall systematic overestimation of bed-elevation (i.e. underestimation of depth) by Bed2013 (see Sect 4.1). The presence of overdeepenings within 25 glacial environments is well established (c.f. Cook and Swift, 2012), their distribution having been observed from bed DEMs for beneath contemporary ice sheets including Greenland (Patton et al., 2016). However, there remains limited quantitative data on their morphology with which to understand the processes of their development (Patton et al., 2015) and the specific relationship between fjord network structure and the locations of overdeepenings and the sills between them. Should additional information become available, such an approach to establish their location could be implemented by introducing rules -for 30 example, an overdeepening of a given step lowering occurs where two fjords of a given width and known depth confluence.
Another approach would be to develop a set of rules which incorporate a fjord hierarchy akin to stream order and their associated Strahler numbers (Strahler, 1957).
The majority of end users of a new Greenland bed DEM including improved bathymetry are expected to be within the ice sheet and polar ocean modelling communities. With this in mind, the approach presented here has been tailored to best suit the purpose of end products that have fjord bathymetry constrained by the synthetic algorithm. Since the algorithm performs better closer to the glacier termini, as opposed to the fjord mouth, users of DEM products based upon this algorithm would be encouraged to focus on processes from the glacier-to-fjord direction (e.g. calving) as opposed to processes focused from the 5 fjord-to-glacier direction (e.g. ocean forcing as in Murray et al., 2010). The impact of high and low frequency stochastic perturbations for topographic datasets for ice sheet modelling is well documented, with models being more sensitive to spatially broad low frequency noise as opposed to higher frequency noise of the same magnitude (Sun et al., 2014). To predict the precise geographical location of sills and overdeepenings with the limited information known for many fjords is a near impossible task.
However, as described in the previous paragraph, the statistical features of these features could be represented by a stochastic 10 model. To the best of our knowledge, our study is the first to consider the statistical properties of fjord bathymetry. This is a significant development as constraining models of high frequency are important where bathymetric surfaces are used to mimic calving (e.g. Lee et al., 2015) or spinning-up ice sheet models over larger regions (e.g. Bindschadler et al., 2013). The exponent for the high frequency tail of the fjord bathymetry power spectrum, 1.74, is consistent with other exponent values found for seafloor topography (Bell, 1975;Goff and Jordan, 1988) and serves as a preliminary guide for future stochastic models. The 15 transition wavelength (∼ 15 km) for the fjord power spectra is shorter than for abyssal hill features in the mid ocean, where the wavelength value is ∼ 40 km (Bell, 1975).
This paper provides a proof of concept routine for constructing geomorphologically realistic fjord geometry in the absence of observations. Actual implementation of the presented routine for large regions (e.g. the Greenlandic coast) would require 20 manual intervention in so far as (i) identifying a seed elevation at the head of the channel and (ii) defining an end zone (e.g. the fjord mouth).
Step (i) could be achieved by using a nearest neighbour approach to acquire the nearest elevation to a given seed location. A solution to step (ii) could be by using an observation density grid where the end zone is identified as being a location with an observation density greater than a chosen value. In addition to this, the values necessary to prevent the development of closed circuit artefacts would have to be adapted to the width of the fjords for which the method is implemented. 25

Summary
Until now, bed-bathymetry DEMs for coastal regions of Greenland have been limited by sparse observations and simplistic interpolation methods being applied within fjord regions. The presented algorithm for synthetic fjord bathymetry provides a new approach to generate bathymetric geometry along fjords. The method takes advantage of observational data where available and assumes that fjords maintain a parabolic cross-sectional profile, thus capturing a leading-order geometric constraint from the 30 ice flow geomorphological processes largely responsible for fjord formation. The validity of the algorithm was tested through comparison with new observational bathymetry data for a fjord system in North West Greenland, and better overall agreement with the data was observed than withfor Bed2013. Additionally, we performed an initial assessment for how a stochastic model of fjord bathymetry could be parameterised, and thus how high frequency perturbations to the flat synthetic geometry could be modelled. The physical validity of the algorithm is limited at multiple channel confluences, as the hierarchy of processes responsible for the landscape features is not explicitly incorporated in the algorithm.
Until more observational data are available, this algorithm provides a suitable estimate for simulating previously unmapped fjord geometry. The presented method will be used to assist in the mapping of fjords within the next Greenland bed DEM data 5 product and has potential application for Antarctica. Using the results of the stochastic model analysis, multiple Greenland bed DEM realisations will be produced, offering the opportunity for the running of ensemble ice sheet model simulations. The release of this new dataset is proposed for 2017.

Data Availability
All data used for the preparation of this manuscript are openly available. The GIMP land classification mask is available and 10 fully documented in Howat et al. (2014). Bed2013 is available and fully documented in Bamber et al. (2013). The IBCAO (v3) DEM is available and fully documented in Jakobsson et al. (2012). The OBS1516 dataset was compiled by I. Fenty and was constructed from (1) the OMG and (2)   Murray, T., Scharrer, K., James, T. D., Dye, S. R., Hanna, E., Booth, A. D., Selmes, N., Luckman, A., Hughes, A. L. C., Cook, S., and Table 1 Parameters and values in the centreline mapping algorithm as applied to the area of interest presented in Fig. 3 Parameter Description Value ∆l finite path segment length 0.8 km r C minimum radius of curvature 6∆l/π ∆a finite path angle difference π/24 l bin bin spatial extent 16 km a bin bin angle extent π/8        Red illustrates where the subtrahend is deeper than the minuend -in the case of OBS156 minus a surface, this indicates an elevation overestimation (too deep).Positive differences (red) occur where the subtrahend is deeper than the minuend, with negative differences (blue) occurring where the subtrahend is shallower than the minuend. Figure 12. Centreline elevation profiles from Bed2013, OBS1516 and the SynthIBCAO and SynthOBS synthetic algorithm approaches. All profiles extend from the head of each fjord to the mouth as depicted in Fig. 3. Figure 13. Composite power spectral density from fjord bed elevation profiles. Blue crosses indicate spectral data that has been averaged over nine fjord profiles, and the solid red curve is the best fit to the parametric model, equation 3. Figure 14. Two different realisations of the stochastic model for high frequency perturbations to the synthetic fjord elevation profiles. The model uses the parametric fit in Fig. 13 to generate the profiles, and is statistically consistent with the OBS1516 bathymetric profiles (green lines in Fig. 12). The overall trend of the fjord bathymetry and lower frequency oscillations (corresponding to wavelengths 14 km or greater) are not synthetically generated and explains why the amplitude of the modelled elevation) is significantly less than the bathymetric observations in Fig 12.