Edinburgh Research Explorer Incorporating modelled subglacial hydrology into inversions for basal drag

. A key challenge in modelling coupled ice-ﬂow– subglacial hydrology is initializing the state and parameters of the system. We address this problem by presenting a workﬂow for initializing these values at the start of a summer melt season. The workﬂow depends on running a subglacial hydrology model for the winter season, when the system is not forced by meltwater inputs, and ice velocities can be assumed constant. Key parameters of the winter run of the subglacial hydrology model are determined from an initial inversion for basal drag using a linear sliding law. The state of the subglacial hydrology model at the end of winter is incorporated into an inversion of basal drag using a non-linear sliding law which is a function of water pressure. We demonstrate this procedure in the Russell Glacier area and compare the output of the linear sliding law with two non-linear sliding laws. Additionally, we compare the modelled winter hydrological state to radar observations and ﬁnd that it is in line with summer rather than winter observations.


Introduction
Subglacial hydrology is an important control on ice velocities at the margin of the Greenland Ice Sheet. Observed seasonal acceleration of ice flow (Joughin et al., 2008;van de Wal et al., 2008;Zwally et al., 2002) is driven by the evolution of the subglacial system between distributed and channelized states in response to meltwater input (Bartholomew et al., 2010;Chandler et al., 2013;Cowton et al., 2013;Schoof, 2010). However, the impact of melt season intensity on seasonal and annual velocities, and how it may change in the future, is not fully understood. Observational studies of land-terminating sectors of the Greenland Ice Sheet reveal a complex set of possible interactions. Increased surface melt may result in faster flow early in the melt season, offset by a stronger late summer deceleration (Sundal et al., 2011). Increased runoff may also lead to more extensive drainage of the ice-sheet base, reducing annual velocities due to slower winter flow . Long-term observations in the ablation zone show surface melt and ice velocities are anticorrelated over decadal timescales (Stevens et al., 2016;Tedstone et al., 2015;van de Wal et al., 2015). The possible impact of surface melt on ice velocities at higher elevations is less well understood, as is the impact in marine-terminating sectors.
Recent subglacial hydrology models have progressed to simultaneously incorporating both distributed and efficient systems, explicitly treating the interaction between the two (de Fleurian et al., 2016;Hewitt, 2013;Schoof, 2010;Werder et al., 2013). These models have shown success in recreating the broad pattern of subglacial development in the summer melt season inferred from GPS measurements (Bartholomew et al., 2011;van de Wal et al., 2015) and dye-tracing experiments Cowton et al., 2013). The development of the subglacial hydrological system has been shown to depend on feedbacks from ice velocities (Hoffman and Price, 2014). However, applications of recent hydrology models coupled with iceflow models have been limited to idealized domains (Hewitt, 2013;Hoffman and Price, 2014;Hoffman et al., 2016;.
Initializing model parameters and state is necessary for applying a linked hydrology-ice dynamics model to the Greenland Ice Sheet. In contrast to the availability of measure-ments at the surface of ice sheets, however, the conditions at the ice-bed interface are poorly constrained. Some key challenges for modelling are the form of the sliding law which relates water pressures to basal drag, the values of the parameters in that relationship, and the values of water pressures at the ice-bed interface.
Inverse methods are an approach which can be used to constrain unknown variables or parameters in an ice-sheet model. Inversions optimize the value of an unknown to minimize the discrepancy between model output and observed data. Since basal drag and the parameters of the sliding law are some of the least constrained inputs to ice-sheet models, a common application of inversions in glaciology is to determine the field of basal drag which best reproduces observed surface velocities. A variety of inversion methodologies have been applied in glaciology. These include iterative methods (Arthern et al., 2015), automatic differentiation (AD; Goldberg and Heimbach, 2013;Heimbach and Bugnion, 2009;Martin and Monnier, 2014), and Lagrangian multiplier methods based on control theory (MacAyeal, 1993;Morlighem et al., 2013).
In this study we develop an ice-sheet model and inversion code, which we apply to the Russell Glacier region of Western Greenland in order to invert for basal drag at the end of winter. The ice-sheet model uses the hybrid formulation of Arthern et al. (2015) and Goldberg (2011) and is numerically similar to Arthern et al. (2015). The inversion procedure is based on AD (Goldberg and Heimbach, 2013). Our application is novel in that we constrain the inversions using the output of a current subglacial hydrology model (Hewitt, 2013), and we investigate the impact of three different sliding laws (linear; Budd as in Budd et al., 1979; and Schoof-type as in Gagliardini et al., 2007, andSchoof, 2005) on the patterns of basal drag predicted by the model. These developments form an important preliminary step towards a fully coupled model for ice dynamics and glacier hydrology which can be validated using current observational ice velocity data and subsequently used for prognostic studies of possible future ice-sheet responses to increased surface melt. The ice-sheet model implemented is based on the hybrid formulation described in Goldberg (2011) and Arthern et al. (2015) and uses the numerical implementation of Arthern et al. (2015). This formulation is derived from the Stokes equations using variational principles (Dukowicz et al., 2010;Goldberg, 2011) and is a hybrid of the shallow ice approximation (Cuffey and Paterson, 2010) and the shallow shelf approximation (MacAyeal, 1989;Morland, 1987). Following Goldberg (2011) and Arthern et al. (2015), the conservation of momentum equations for depth-averaged velocities are where u(x, y, z) and v(x, y, z) are velocities in the x and y directions, η(x, y, z) is dynamic viscosity, h(x, y) is ice thickness, s(x, y) is surface elevation, τ bx(x,y) and τ by (x, y) are basal drag in the x and y directions, g is the magnitude of gravitational acceleration, and ρ i is the density of ice. The overbar (x) denotes the depth-averaged value of a variable, so thatū(x, y) andv(x, y) are depth-averaged velocities and η(x, y) is depth-averaged viscosity. Basal drag is defined by the sliding law. Three different sliding laws are implemented in the ice-sheet model: where τ b = (τ bx (x, y), τ by (x, y)) is the basal drag, u(x, y, b), v(x, y, b)) is the basal velocity, U b is the sliding speed (|u b |), N (x, y) = ρ i gh−p w is the effective pressure at the ice-sheet bed, p w is water pressure, β(x, y) is a basal drag coefficient, µ a (x, y) is a drag coefficient, p and q are positive exponents, µ b (x, y) is a limiting roughness slope, λ b is the characteristic bed roughness length, and A b and n are coefficients in Glen's flow law (Hewitt, 2013). A b is the ice creep parameter used for basal ice. It set an order of magnitude lower than A to account for warmer ice at the base (following Hewitt, 2013). The value of A b is used in both the Schoof sliding law and the subglacial hydrology model for determining creep closure of channels and cavities. The value of A is used in the momentum equations. Following Hewitt (2013), negative effective pressures are eliminated by setting N + = max(N, 0) and regularized with a small regularization constant.
The linear sliding law (Eq. 3) represents all ice-bed interactions by a single friction coefficient β. The second and third equations are a Budd sliding law and a Schoof sliding law respectively. These attempt to explicitly represent more complex interactions at the ice-bed interface, in particular the impact of basal water pressure. Equation (4) is a power law commonly used in glaciology to describe basal rheology (e.g. Bueler and Brown, 2009;MacAyeal, 1989;Hewitt, 2013), although typically with no dependence on effective pressure (p = 0). At high effective pressures the Schoof sliding law has a similar form ( It is useful to represent the sliding laws in a common form: where C is a function multiplying basal velocities. The form and parameters of C depend on the sliding law.
The boundary conditions at the terminating margin of the ice sheet are 2ηh(2∂ yv + ∂ xū )n y +ηh(∂ yū + ∂ xv )n where ρ w is the density of water, d is the ice draft (zero at land-terminating portions of the margin), andn x andn y are the components of the outward pointing unit vector normal to the terminating margin (Goldberg, 2011;Arthern et al., 2015). Two further boundary conditions are used in the ice-sheet model: a no-penetration condition at the margin of nunataks and a Dirichlet boundary condition at the lateral margins of the ice-sheet domain which are not the termination edge.
The equation for viscosity is where 0 is a regularization term. Vertical shearing in the hybrid formulation is approximated by As in Goldberg (2011) and Arthern et al. (2015), a linear relationship between vertical shear stresses and depth is assumed: Viscosity is defined implicitly by Eq. (9). With the standard choice of n = 3, this is a cubic equation and can be solved exactly. Alternatively, a previous value of viscosity can be used to calculate an updated value. This process can be iterated upon to create a fixed point iteration. The default procedure in the model is to do two iterations (Koziol, 2017).
The hybrid formulation of the conservation of momentum equations depend on depth-integrated viscosity: This integral, and others, is numerically integrated using the composite Simpson law.
Following Arthern et al. (2015), the following integral is defined: This integral can be used to define expressions for surface velocity in terms of basal velocity and basal velocity in terms of depth-averaged velocity (Arthern et al., 2015): where F 1 and F 2 are determined using Eq. 13 .
Additionally, defining C eff as follows, leads to an expression for basal drag in terms of depthaveraged velocity (Goldberg, 2011;Arthern et al., 2015):

Model implementation
As in Arthern et al. (2015), Eqs. 1 and 2 can be written in the following form: where and f = ρ i gh∂ x s ρ i gh∂ y s .
Equation (18) is a non-linear equation for depth-integrated velocity. The non-linearity arises since depth-integrated viscosity is a function of velocity and, in the case of a non-linear www.the-cryosphere.net/11/2783/2017/ The Cryosphere, 11, 2783-2797, 2017 sliding law, since C eff is also a function of velocity. The icesheet model solves Eq. 18 on an Arakawa-C finite difference grid using a Picard iterative process. Equation (18) is discretized following Arthern et al. (2015). The primary difference is that operators are appropriately extended to allow for periodic boundary conditions in the ISMIP-HOM experiments (Pattyn et al., 2008). Discretization of Eq. 18 results in a linear system of equations, which can be written as where the matrix (L) corresponds to the operator L, while the vector x corresponds toū, and the vector b corresponds to f . MATLAB's backslash operator is used to solve this system of equations. Alternatively, preconditioned iterative methods can be used (Arthern et al., 2015;Goldberg and Heimbach, 2013). The Picard iteration linearizes Eq. 18 by constructing L using the velocity of the previous iteration. An initial velocity guess and viscosity guess form the initial L. Equation (18) is then solved for an updated velocity guess, which in turn can is used to update viscosity and C eff . This process is repeated within a loop until the solution converges below a specified tolerance or until a prescribed number of iterations are reached.
Evolution of surface geometry is not included in the icesheet model. This is appropriate since the ice-sheet model is applied on annual timescales, over which significant changes in ice-sheet geometry are not expected.
The ice-sheet model was tested against the ISMIP-HOM benchmark experiments A and C (Pattyn et al., 2008) and found to compare favourably against previous models (Koziol, 2017).

Model formulation
This section describes the details of an inversion code developed in conjunction with the ice-sheet model. The methodology is based on Goldberg and Heimbach (2013). However, the implementation developed here has a more limited capability due to software limitations.
The cost function returns a scalar which measures the fit of the model to the observations. The cost function is defined as where γ 1 and γ 1 are user-defined scaling factors, s is the surface domain, b is the basal domain, w(x, y) is a weighting function, U obs (x, y) are observed surface ice speeds, U s (x, y) are modelled surface speeds, and α(x, y) is the control parameter.
The cost function defined above has two terms: J = γ 1 J 0 + γ 2 J Reg . The first term (J 0 ) measures the weighted square of the difference between observed and modelled velocity. The second term (J Reg ) is a Tikhonov regularization term, which penalizes oscillations in α and stabilizes the inversion (Morlighem et al., 2013). Other formulations of the cost function are possible (e.g. Morlighem et al., 2013).
The control parameter refers to the variable which the inversion process optimizes in order to best match model prediction and observations. Since our aim is to determine the basal drag, the control parameter is a parameter in the basal sliding law. For the linear sliding law, α = β 2 . For the Budd sliding law, α = µ a . Although the Schoof sliding law has two unknowns which can be inverted for, µ b exerts a dominating control. Hence, λ b is set to a constant while α = µ b . In the numerical implementation of the adjoint, α is parameterized as α(x, y) = exp ζ (x, y) . This ensures that α remains positive, as expected for each of the three sliding laws. For simplicity, this is neglected in the remainder of the paper, and the discussion focuses on recovering α rather than ζ .
The inversion process aims to determine the field of α which minimizes the cost function. This is an optimization problem. Starting with an initial guess for α, the gradient of the cost function with respect to the initial input α is determined. The gradient provides a search direction for the optimization algorithm, which updates α. This process is repeated iteratively until α converges below a tolerance or until a maximum number of iterations occur. The critical component in this process is the gradient dJ dα . The process to calculate this gradient is described in the next two subsections.

Adjoint model description
The methodology to obtain the gradient dJ dα follows from Goldberg and Heimbach (2013). The key concepts of this approach are first explained for a generic algorithm, before showing how they can be applied to the ice-sheet model. This explanation follows that of Errico (1997) and Goldberg and Heimbach (2013). Consider the following model: where φ is an arbitrary variable (or array of variables), and B can be considered a sequence of operations: and each operation can be written as b N = B N (). Further, define a function J: where J returns a scalar. In the context of the adjoint model, the function is known as the cost function, objective function, or target function (Goldberg and Heimbach, 2013). This function quantifies an aspect of the model output which is of interest, such as the mean error of model output relative to observations. The aim is to determine the gradient of the cost function J with respect to the initial input φ. To provide context for the adjoint model, the tangent linear model (TLM) is presented first. In the TLM, a small perturbation in the input is propagated forward through the model to determine the corresponding perturbation in the output. Applying the chain rule to J = J (b) = J (B(φ)) leads to the corresponding TLM: There are several observations about the TLM. First, the TLM determines the perturbation of δJ from the perturbation of a single element φ i . As the perturbation δφ i approaches zero, δJ δφ i converges to dJ dφ i . Second, to determine dJ dφ , the TLM needs to be run for each entry in φ. Although for small models this approach is feasible, the computational cost is too great for glaciological problems on domains of the size of interest. Finally, the TLM acts in a similar direction as the model B, in that the functions are applied successively starting with the counterpart to B 0 (Errico, 1997).
The concept behind the adjoint model is that rather than determining how changes in the input φ impact the cost function J , it can be more efficient to determine how changes in the cost function J impact the initial input φ. In the adjoint model, sensitivities of J are propagated backwards through the model to determine the resulting change in φ. Similar to the TLM, the adjoint model is derived by applying the chain rule to J = J (b) = J (B(φ)): Key observations about the adjoint model are as follows.
(1) In contrast to the TLM, which acts upon a perturbation, the adjoint model acts upon the sensitivity of the cost function.
(2) A single run of the adjoint model is sufficient to determine the gradient δJ δφ . (3) The adjoint model runs in reverse relative to both the model and the TLM, in that the adjoint model applies functions beginning with the counterpart to B N and ending with the counterpart of B 0 (Errico, 1997).

Adjoint model implementation
The adjoint model is generated based on AD (Griewank and Walther, 2008) of the MATLAB code implementations of the forward model. AD tools process an input code to generate a counterpart code, which returns the corresponding gradient (or Jacobian). The central concept behind AD is that a computer program is fundamentally a sequence of elementary operations and functions. This admits the repeated application of the chain rule to generate a derivative of high accuracy.
Multiple methodologies exist for AD tools to generate the derivative code. Previous applications of AD software to generate the adjoint in glaciology (Heimbach and Bugnion, 2009;Goldberg and Heimbach, 2013;Martin and Monnier, 2014) have used reverse accumulation AD tools (e.g. Giering et al., 2005;Hascoet and Pascual, 2004). These types of AD software are conceptually similar to the adjoint model. They are designed to determine the gradient of a function (input code) by propagating sensitivities of the output variables backwards to the input variables. Hence, an ice-sheet model can be processed with relatively little modification by reverse accumulation AD tools to generate the adjoint model.
Here, we apply the open source AD tool ADiGator (Weinstein andRao, 2011-2016), which in contrast to previous work is a forward accumulation AD tool. The methodology of forward accumulation is conceptually similar to the TLM. It is designed to determine the gradient of a function (input code) by propagating sensitivities of the input variables forward through the program to the output. Applying a forward AD tool to an ice-sheet model to generate the adjoint is not feasible due to the size of the control space. Rather, we generate the adjoint by applying ADiGator to segments of the ice-sheet model code and multiplying the resulting Jacobians following Eq. (27).
Pseudocode of the main ice-sheet model routine is shown in Algorithm A1, and the corresponding code to calculate the adjoint is shown in Algorithm A2 (see Appendix). Two new functions, S1 and S2, appear in the adjoint code. These encapsulate segments of code from the forward model and can be processed by ADiGator. The function S2 contains code which spans over two Picard iterations. The adjoint does not contain a for loop corresponding to iterating through the Picard iterations in reverse (cf. Goldberg et al., 2016). Rather, values from the final two Picard iterations of the forward model are saved and used as input for the adjoint code. The adjoint model is also modified to solve the cubic equation (following Arthern et al., 2015) to determine η, rather than storing values from the previous iterations and implementing a fixed point iteration. This impacts the η,η, and F a functions but leaves the overall structure the same. This is a necessary modification for ADiGator.
The adjoint code explicitly calculates several Jacobian matrices (lines 15 to 23 in Algorithm A2). ADiGator is applied to the corresponding functions to generate the Jacobian matrices, except the solution to the system of linear equations, which requires special treatment. A counterpart to the linear solve which returns the corresponding derivate is manually programmed following the procedure detailed in the appendix of Martin and Monnier (2014). The adjoint is then calculated by multiplying out the sensitivities of the cost function with the transposes of the Jacobian matrices. Although this process is more complicated and less flexible than previous approaches, it is necessary because no noncommercial AD reverse accumulation tool is available for MATLAB. This implementation of the adjoint is equivalent to previously published adjoint implementations (Goldberg and Heimbach, 2013;Martin and Monnier, 2014) restricted to one reverse step in the Picard iteration. This is mathematically equivalent to the Lagrangian multiplier method introduced by MacAyeal (1993) (Heimbach and Bugnion, 2009).
The gradient from the adjoint model is used to solve the optimization problem which minimizes the cost function. The inversion code relies on minFunc (Schmidt, 2005), a publicly available MATLAB unconstrained optimization package. The L-BFGS routine, with a Wolfe condition backtracking line search, is applied in the inversion code. The cost function is discretized using the same finite difference operators as the ice-sheet model.
Performance of the inversion code was verified using a series of identical twin tests (Goldberg and Heimbach, 2013). Results are shown in Koziol (2017).

Subglacial hydrology model
This subglacial hydrology model used is described in detail in Hewitt (2013) and Banwell et al. (2016) and is similar conceptually to the model presented in Werder et al. (2013). Here, the version employed in Banwell et al. (2016) is applied.
Both distributed and channelized flow are represented in the subglacial hydrology model. Distributed flow is described by an average thickness and flux over a representative area. As in Banwell et al. (2016), the distributed system is composed of two components: a cavity sheet and an elastic sheet. The elastic sheet is included to simulate "hydraulic jacking" from lake hydrofracture events and is activated only when the effective pressure drops to zero and below. Channels have the potential to form along the edges and diagonals of the numerical grid. Channels are initiated by dissipative heating from the distributed system over an incipient channel width length scale. The model is written in MATLAB, using a finite difference numerical grid and an implicit forward time step method. For full details, consult Hewitt (2013) and Banwell et al. (2016).

Application to Russell Glacier area
The Russell Glacier area is a land-terminating sector of the Greenland Ice Sheet (Fig. 1). The ice-sheet model and inversion code are applied to the Russell Glacier area to determine the basal boundary condition at the end of the winter season.
An outline of the study area is shown in (Fig. 1). The northern and southern boundaries are selected to be roughly in line with basal watersheds determined using the Shreve (1972) approximation for hydraulic gradient. The northern boundary is approximately the same as used by Bougamont et al. (2014) andde Fleurian et al. (2016). The southern boundary is further south relative to Bougamont et al. (2014) but north of the southern boundary in de Fleurian et al. (2016). The eastern boundary was selected to extend up ice of the GPS stations (Tedstone and Neinow, 2017). The western boundary is the ice margin. There is a nunatak near the western boundary.
The ice-sheet model-inversion code is applied to determine the basal boundary condition at the end of the 2008-2009 winter season in the Russell Glacier study site. The end of the winter season is assumed to be day 120 of the year (30 April). Although the exact day is somewhat arbitrary, this day was selected as it is shortly before surface runoff begins in the study area and shortly before GPS records in the study site show enhanced motion (Tedstone and Neinow, 2017). Applying the ice-sheet model-inversion code to the Russell Glacier area requires a number of datasets. Mean winter surface velocities for 2008/2009 (Fig. 2) are provided by the MEaSUREs Greenland Ice Sheet Velocity Map at 500 m resolution (Joughin et al., 2010a, b). Surface and basal topography (Fig. 2) is provided by the BedMachine2 dataset (Morlighem et al., 2014(Morlighem et al., , 2015 and is resampled to 500 m resolution from 150 m resolution to match the velocity data. This is slightly coarser than the reported true resolution of 400 m for the ice thickness. The 500 m grid resolution results in a grid size of 132 × 274 for the domain. Fifty vertical layers are used for integration using Simpson's rule. An important assumption made is that the mean winter velocities are representative of both the beginning and end of winter. This assumption is justified by observing published GPS records in southwest Greenland (Colgan et al., 2012;van de Wal et al., 2015). These observations show that although velocities increase throughout the winter, the magnitude of the change is limited (<25%).
Inversions are initialized using a basal drag set to the local driving stress smoothed by a 3 × 3 grid cell mean filter. The ice-margin boundary is described in the ice-sheet model by Eqs. (7) and (8) while on the three other boundaries a Dirichlet boundary condition is applied. The inverse values of the errors provided with the surface velocity measurements are used as weights in the cost function.
The results of inversions depends on the relative values of the scaling factors γ 1 to γ 2 in the cost function (Eq. 22). For each sliding law, a series of inversions is performed with γ 1 set to 1 while varying γ 2 . An L-curve analysis is applied to select the inversion which best balances fitting the velocity observations while penalizing spurious oscillations in basal drag. Parameters for the ice-sheet model-inversion code are listed in Table 1. Similar to Hewitt (2013), the ice-flow creep parameter (A) is selected to be 7·10 −25 Pa −3 s −1 . This corresponds to an ice temperature of approximately −7 • C (Cuffey and Paterson, 2010). This choice for A results in the ratio of basal velocity to surface velocity remaining greater than 0.5 throughout the study area.
The parameters for the subglacial hydrology are the result of an extensive parameter search using a coupled iceflow-subglacial hydrology model in Koziol (2017). Many of the parameters are the same as published in Banwell et al. (2016) andHewitt (2013). However, testing of the reported optimal parameters for the Paakitsoq region reported by Banwell et al. (2016) using the integrated model showed poor agreement with GPS measurements due to insufficient volumes of water being evacuated from mid-elevations.
A workflow is developed for incorporating modelled effective pressure into inversions using non-linear sliding laws. This workflow is motivated by the idea that both the sub-  glacial hydrological system and ice flow are in quasi-steady state during the winter. This allows us to invert for background values of the constants in the sliding laws. The initial step is to invert using a linear sliding law for the basal drag coefficient. Basal velocities are calculated from modelled depth-integrated velocities (Eq. 15). The modelled basal drag and basal velocities then provide the necessary input for the subglacial hydrology model to calculate a distributed basal melt rate. The modelled distributed basal melt rate incorporates geothermal flux but neglects heat loss to the interior of the ice sheet (Hewitt, 2013). The subglacial hydrology model is then run for the winter season with the basal drag and basal velocities from the linear inversion. The model is run at 500 m resolution (identical to the inversions), with no-flow boundary conditions at the northern, southern, and eastern boundaries. The ice margin is assumed to be at atmospheric pressure. This boundary condition is modified at necessary places to prevent inflow of water from beyond the ice-sheet margin. Similarly to Banwell et al. (2016), the subglacial hydrology model is initialized with the thickness of the sheet flow layer set to 0.10 m. Testing showed that varying initial thickness has negligible impact. At this stage, the ice-sheet model remains unconnected, and the input basal velocities are assumed to be constant. The subglacial hydrology model run provides a modelled water pressure distribution over the study site.
Finally, the non-linear inversions are run using the modelled water pressure from the subglacial hydrology model winter run. Two sets of inversions are conducted, one for the Budd sliding law and one for the Schoof sliding law. The first set of inversions seeks to determine the distribution of µ a , while the second inverts for µ b . Similar to the linear sliding law, an L-curve analysis is employed to determine the relative values of γ 1 to γ 2 .

Linear inversion
Six inversions using the linear sliding law are run (Fig. 3). Using the L-curve plot, the inversion with γ 2 = 1 · 10 −12 is selected as optimal. The value of J 0 for this inversion is 1.56· 10 11 .
A map of the difference between observed and modelled velocities shows the highest difference occurs along the ice margin and in the vicinity of the nunatak (Fig. 4). Figure 5 shows the inverted basal drag parameter, basal drag, and the sliding ratio ( U b U s ) for the linear sliding law.

Subglacial hydrology model
Basal melt during the winter is shown in Fig. 6. Most values are between 0.015 and 0.03 m yr −1 , with higher values predominately occurring near the nunatak. The spatial pattern of melt broadly reflects the patterns of surface velocities (Fig. 2). The subglacial hydrology model winter run evolves rapidly at the beginning of the run. However, by day 50 of the model run the rate of change is significantly reduced, and by day 240 of the run the model is in an approximate steady state.
The distribution of sheet thickness at the end of winter mirrors basal topography, with the sheet thickest in topographic lows (Fig. 7). The maximum distributed system sheet thickness is 0.36 m. The effective pressure also reflects the basal topography, with lowest effective pressures located in topo-graphic lows. Since the lowest effective pressure is 0.44 MPa, no part of the ice sheet is near flotation. The model predicts minor channelization in two locations (not shown), with single channels extending from the margin several kilometres.

Non-linear sliding laws
An L-curve analysis (Fig. 3) is used to determine the optimum inversion for each of the non-linear sliding laws. The inversions corresponding to γ 2 = 1 are selected for the Budd sliding law, while the inversion corresponding to γ 2 = 10 11 is selected for the Schoof sliding law. These were selected so that the cost term of the inversions were similar to that of the linear sliding law. The two cost terms for the Budd and Schoof sliding laws are J 0 = 1.78 · 10 11 and J 0 = 1.60 · 10 11 respectively.   Figure 5 shows the inversion results from the Budd and Schoof sliding laws. Inverted basal drag and the consequent sliding ratio using the two non-linear sliding laws are very similar to the results from the linear sliding law. Model mismatch is also similar for all three sliding laws (Fig. 4).

Discussion
Inversions of the Russell Glacier area are run with a constant creep parameter A, corresponding to an ice temperature of approximately −7 • C (Cuffey and Paterson, 2010). For inversions with the linear sliding law, tests showed poorer results when A increased (corresponding to warmer ice). As A decreased, the sliding ratio approached one uniformly, and inversion results were better able to match observed surface velocities. The value of A was selected as a balance of model fit, while keeping a contribution to motion from internal de-formation. Observations from two boreholes located in the Paakitsoq region show that internal deformation results in approximately 27-56 % of ice velocities during winter (Ryser et al., 2014). In reality, A would have a heterogeneous distribution. By using a constant A, the basal drag parameter will account for some of the effects, which would otherwise be due to variation in A.
Basal velocities determined from the optimal inversion using a linear sliding law are input into the subglacial hydrology model. The distribution of basal velocities is used to calculate both the basal melt rate and the cavity space in the continuum sheet flow. Due to the selection of a creep parameter such that the sliding ratio is relatively high, it is likely that basal velocities are overestimated. This would result in an overestimate of water generated at the ice-bed interface and an overestimate of the capacity of cavity space. Application of a higher order ice-sheet model would be advantageous in these regards.
The pattern of basal drag inverted using the three different sliding laws show limited differences. This is due to the fact that basal shear traction must satisfy the global stress balance (Joughin et al., 2004;Minchew et al., 2016). The basal drag and basal velocities from the linear sliding law used to initiate the subglacial hydrology model are therefore self consistent with the subsequent inversion results of the two non-linear sliding laws. For the winter effective pressures predicted by the subglacial hydrology model, we find that the Schoof sliding law is in the viscous drag regime. For a representative basal velocity of 75 m yr −1 , the transition to Coulomb friction occurs at effective pressures of approximately 0.7 MPa. This is below modelled effective pressures, which are above 1.3 MPa for most of the study domain.
Interpretation of radar lines in the Russell Glacier area suggests significant winter storage of water along topographic highs, while significant water flow through topographic lows occurs during the summer melt seasons (Chu et al., 2016). Based on these observations, the subglacial hydrology runs are reflective of summer conditions rather than winter conditions. Water storage, which would be characterized by high sheet thickness, is not observed along topographic highs. Chu et al. (2016) attribute storage on topographic ridges to water storage in parts of the distributed system which become isolated at the end of the melt season. In contrast, porous sediments in bedrock troughs are hypothesized to allow water to drain (Chu et al., 2016). The treatment of the bed in the subglacial hydrology model is uniform. It does not account for differences in till cover or bed properties, nor does it account for sub-grid-scale heterogeneity in the distributed system, which is likely the cause of water storage. Replicating these observations likely requires the implementation of another model component, such as the weakly connected distributed system proposed by Hoffman et al. (2016). In general, model output from the subglacial hydrology model can be expected to be much more sensitive to the model formulation during the winter than the summer, when the system is forced by high water input. In line with inferences from tracer injections , the model does not predict a channelized system at the margin during the winter.
The initialization procedure introduced is not capable of producing the inferred year-on-year differences in the subglacial hydrological system at the end of winter. Currently, the subglacial hydrology reaches an approximate steady state by day 240 and is not particularly sensitive to the initialization of the distributed sheet thickness. A full steady state takes approximately 2 years (Hewitt, 2013). In contrast, observations suggest that summer melt has an impact on the state of the hydrological system during the subsequent winter (Chu et al., 2016;Sole et al., 2013). The model output therefore can only be considered an approximation to a generic hydrological state. Any discrepancy between the modelled and actual hydrological system is expected to have a greater impact on inversions using the Schoof sliding law, since it has a stronger dependence on effective pressure. In the limit of viscous flow, the Schoof sliding law depends on N. In contrast, the Budd law is a function of N 1/3 (Budd et al., 1979). All inversions are conducted using mean winter velocities from 2008 to 2009. Annual differences in mean winter velocities are expected to have a minimal impact, as observed year-on-year differences are on the order 20 m yr −1 , which is not significantly greater than the velocity mismatch in the inversions.
Other procedures for determining the background parameters of sliding laws can likely be devised. Currently the procedure only uses mean winter velocities. Using mean annual velocities may improve estimates of the sliding law parameters by incorporating information from the melt season. A subglacial hydrological model could be run for an entire year and basal parameters determined from an annual average water pressure. A key difficulty is running the hydrological model during the summer, as the development of the system is known to depend on feedbacks with velocity (Hoffman and Price, 2014). This issue can be avoided by using velocity measurements from remote sensing as a model forcing (e.g Fahnestock et al., 2016). An advantage of running the subglacial hydrology model during the summer months is that model output may be more representative of water flow beneath the ice sheet. Although in its current form the model is too complex, a simplified subglacial hydrology model may be suitable to time-dependent adjoint modelling (Goldberg and Heimbach, 2013). Here, we have assumed that the parameters of the sliding law are independent of time. This assumption is better suited for bedrock than till, as properties of till are dependent on saturation and deformational history (Minchew et al., 2016).

Conclusions
A new ice-sheet model and adjoint code are presented. The ice-sheet model is coupled to a recent subglacial hydrology model (Hewitt, 2013). A procedure for initializing a coupled subglacial hydrology-ice-sheet model using a winter run is also proposed. The modelled state of the subglacial hydrological system at the end of winter appears to reflect summer observations rather than winter observations. This is likely the result of model formulation rather than the initialization procedure and could lend support to the potential need for an additional, weakly connected component in hydrological models (Hoffman et al., 2016). However, the initialization procedure presented here will continue to prove useful as model development advances, as this is independent of the hydrological model used. The results are subsequently used to run inversions using non-linear sliding laws which are functions of effective pressure. This allows the background parameters for the sliding law to be determined. To date, this appears to be the first work to incorporate modelled water pressures in an inversion and the first to invert with a sliding law explicitly dependent on effective pressure. The usefulness of this inversion for initiating coupled ice-sheethydrology model simulations is shown in an upcoming publication.
Code and data availability. All datasets used are publicly available. BedMachine data are available from Morlighem et al. (2015). MEaSUREs Greenland Ice Sheet Velocity Map data are available from Joughin et al. (2010b). Code is currently not available.