Journal cover
Journal topic
**The Cryosphere**
An interactive open-access journal of the European Geosciences Union

Journal topic

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications
- Imprint
- Data protection

- About
- Editorial board
- Articles
- Special issues
- Highlight articles
- Manuscript tracking
- Subscribe to alerts
- Peer review
- For authors
- For reviewers
- EGU publications
- Imprint
- Data protection

**Research article**
07 Mar 2018

**Research article** | 07 Mar 2018

Simulating the thermomechanical features of Laohugou Glacier No. 12

^{1}State Key Laboratory of Cryospheric Science, Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou 730000, China^{2}University of Chinese Academy of Sciences, Beijing 100049, China^{3}Institute of Polar Meteorology, Chinese Academy of Meteorological Sciences, Beijing 100081, China^{4}College of Geography and Environment, Shandong Normal University, Jinan 250014, China^{a}now at: Fluid Dynamics and Solid Mechanics Group, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

Abstract

Back to toptop
By combining in situ measurements and a two-dimensional thermomechanically coupled ice flow model, we investigate the thermomechanical features of the largest valley glacier (Laohugou Glacier No. 12; LHG12) on Qilian Shan located in the arid region of western China. Our model results suggest that LHG12, previously considered as fully cold, is probably polythermal, with a lower temperate ice layer overlain by an upper layer of cold ice over a large region of the ablation area. Modelled ice surface velocities match well with the in situ observations in the east branch (main branch) but clearly underestimate those near the glacier terminus, possibly because the convergent flow is ignored and the basal sliding beneath the confluence area is underestimated. The modelled ice temperatures are in very good agreement with the in situ measurements from a deep borehole (110 m deep) in the upper ablation area. The model results are sensitive to surface thermal boundary conditions, for example surface air temperature and near-surface ice temperature. In this study, we use a Dirichlet surface thermal condition constrained by 20 m borehole temperatures and annual surface air temperatures. Like many other alpine glaciers, strain heating is important in controlling the englacial thermal structure of LHG12. Our transient simulations indicate that the accumulation zone becomes colder during the last two decades as a response to the elevated equilibrium line altitude and the rising summer air temperatures. We suggest that the extent of accumulation basin (the amount of refreezing latent heat from meltwater) of LHG12 has a considerable impact on the englacial thermal status.

How to cite

Back to top
top
How to cite.

Wang, Y., Zhang, T., Ren, J., Qin, X., Liu, Y., Sun, W., Chen, J., Ding, M., Du, W., and Qin, D.: An investigation of the thermomechanical features of Laohugou Glacier No. 12 on Qilian Shan, western China, using a two-dimensional first-order flow-band ice flow model, The Cryosphere, 12, 851-866, https://doi.org/10.5194/tc-12-851-2018, 2018.

1 Introduction

Back to toptop
The storage of water in glaciers is an important component of the
hydrological cycle at different time scales
(Huss et al., 2010; Jansson et al., 2003), especially in arid and semiarid regions
such as northwestern China, where many glaciers are currently
retreating and disappearing
(Neckel et al., 2014; Tian et al., 2014; Yao et al., 2012). Located on the northeastern edge
of the Tibetan Plateau (36–39^{∘} N, 94–104^{∘} E),
Qilian Shan (QS) develops 2051 glaciers covering an area of
approximately 1057 km^{2} with a total ice volume of
approximately 51 km^{3} (Guo et al., 2014, 2015). Meltwater from
QS glaciers is a very important water resource for the agricultural
irrigation and socioeconomic development of the oasis cities in
northwestern China. Thus, the changes in the QS glaciers that occur
as the climate becomes warmer in the near future are concerning.

Due to logistic difficulties, few QS glaciers have been investigated
in previous decades. However, Laohugou Glacier No. 12 (hereafter
referred to as LHG12), the largest valley glacier of QS, has been
investigated. Comprised of two branches (east and west), LHG12 is
located on the north slope of the western QS
(39^{∘}27^{′} N, 96^{∘}32^{′} E;
Fig. 1), with a length of
approximately 10 km, an area of approximately 20 km^{2},
and an elevation range of 4260–5481 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$
(Liu et al., 2011). LHG12 was first studied by a Chinese expedition
during 1958–1962 and was considered again in short-term field
campaigns in the 1970s and 1980s that were aimed at monitoring glacier
changes (Du et al., 2008). Since 2008, the Chinese Academy of Sciences
has operated a field station for obtaining meteorological and
glaciological measurements of LHG12.

The temperature distribution of a glacier primarily controls the ice rheology, englacial hydrology, and basal sliding conditions (Blatter and Hutter, 1991; Irvine-Fynn et al., 2011; Schäfer et al., 2014). A good understanding of the glacier thermal conditions is important for predicting glacier response to climate change (Gilbert et al., 2015; Wilson et al., 2013), improving glacier hazard analysis (Gilbert et al., 2012), and reconstructing past climate histories (Gilbert et al., 2010; Vincent et al., 2007). Different energy inputs at the glacier surface determine the near-surface ice temperature. In particular, the latent heat released by meltwater refreezing in the percolation zone can largely warm the near-surface area (Blatter, 1987; Blatter and Kappenberger, 1988; Gilbert et al., 2014b; Huang et al., 1982; Müller, 1976). The importance of surface thermal boundary condition in controlling the thermal regime of a glacier has been recognized by recent thermomechanical glacier model studies (Gilbert et al., 2014a; Meierbachtol et al., 2015; Wilson and Flowers, 2013). Using both in situ measurements and numerical models, Meierbachtol et al. (2015) argued that shallow borehole ice temperatures served as better boundary constraints than surface air temperatures in Greenland. However, for the east Rongbuk Glacier on Mt. Everest, which is considered polythermal, Zhang et al. (2013) found that the modelled ice temperatures agreed well with the in situ shallow borehole observations when using surface air temperatures as the surface thermal boundary condition. Therefore, careful investigation of the upper thermal boundary condition is highly necessary for glaciers in different regions under different climate conditions.

LHG12 is widely considered as an extremely continental type (cold) glacier and is characterized by low temperatures and precipitation (Huang, 1990; Shi and Liu, 2000). However, in recent years, we have observed extensive and widespread meltwater on the ice surface and at the glacier terminus. In addition, percolation of snow meltwater consistently occurs in the accumulation basin during the summer. Therefore, we address two pressing questions in this study. Firstly, what is the present thermal status of LHG12? Secondly, how do different surface thermal boundary conditions impact the modelled ice temperature and flow fields? Because temperate ice can assist basal motion and accelerate glacier retreat, understanding the current thermal status of LHG12 is very important for predicting its future dynamic behaviour.

To answer these questions, we conduct diagnostic and transient simulations for LHG12 by using a thermomechanically coupled first-order flow-band ice flow model. This paper is organized as follows: first, we provide a detailed description of the glaciological datasets of LHG12. Then, we briefly review the numerical ice flow model and the surface mass balance model used in this study. Next, we perform both diagnostic and transient simulations. In the diagnostic simulations, we first investigate the sensitivities of ice flow model parameters, and we then explore the impacts of different surface thermal boundary conditions and assess the contributions of heat advection, strain heating, and basal sliding to the temperature field of LHG12. Transient simulations are performed during the period of 1961–2011. We then compare the transient results with measured ice surface velocities and ice temperature profile obtained from a deep borehole. We also investigate the evolution of the temperature profile in the deep borehole and the impacts of initial equilibrium line altitude (ELA) on the thermal regime of the glacier. Finally, we discuss the limitations of our model and present the important conclusions that resulted from this study.

2 Data

Back to toptop
Most in situ observations, e.g. borehole ice temperatures, surface air temperatures and ice surface velocities, have been made on the east branch (main branch) of LHG12 (Fig. 1). Measurements on the west branch are sparse and temporally discontinuous. Thus, we only consider the in situ data from the east tributary when building our numerical ice flow model.

In July–August 2009 and 2014, two ground-penetrating radar (GPR) surveys were conducted on LHG12 using a pulseEKKO PRO system with centre frequencies of 100 MHz (2009) and 50 MHz (2014) (Fig. 1b). Wang et al. (2016) have presented details regarding the GPR data collection and post-processing.

As shown in Fig. 2, the east branch of LHG12 has
a mean ice thickness of approximately 190 m. We observed the
thickest ice layer (approximately 261 m) at
4864 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Generally, the ice surface of LHG12 is gently
undulating, with a mean slope of 0.08^{∘}, and the bed of LHG12
shows significant overdeepening in the middle of the centre flow line
(CL) (Fig. 2). To account for the lateral effects
exerted by glacier valley walls in our two-dimensional (2-D) ice flow model, we
parameterize the lateral drag using the glacier half widths. Based on
the GPR measurements on LHG12, we parameterize the glacier
cross sections by a power law function *z*=*a**W*(*z*)^{b}, where *z* and
*W*(*z*) are the vertical and horizontal distances from the lowest point
of the profile, and *a* and *b* are constants representing the
flatness and steepness of the glacier valley, respectively
(Svensson, 1959). The *b* values for LHG12 range from 0.8 to 1.6,
indicating that the valley is approximately V-shaped
(Wang et al., 2016). As an input for the flow-band ice flow model, the
glacier width, *W*, was also calculated by ignoring all tributaries
(including the west branch) (Figs. 1b and
2).

The ice surface velocities of LHG12 were determined from repeated
surveys of stakes drilled into the ice surface. All stakes were
located in the distance between 0.6 and 7.9 km along the CL
(Fig. 3), spanning an elevation range of
4355–4990 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. 1b). We measured
the stake positions using a real-time kinematic (RTK) fixed solution
by a South Lingrui S82 GPS system (Liu et al., 2011). The accuracy of the
GPS positioning is on the order of a few centimetres and the
uncertainty of the calculated ice surface velocities is estimated to be
less than 1 m a^{−1}. Because it is difficult to conduct
fieldwork on LHG12 (due to e.g. crevasses and supraglacial
streams), it was nearly impossible to measure all stakes each
observational year. Thus, the current dataset includes annual ice
surface velocities from 2008 to 2009 and from 2009 to 2010, summer measurements
from 17 June to 30 August 2008, and winter measurements from
1 February to 28 May 2010.

The in situ ice surface velocities shown in Fig. 3 are all
from stakes near the CL (Fig. 1b). Small ice surface
velocities (<17 m a^{−1}) are clearly visible in the upper
accumulation (0–1.2 km) and lower ablation areas (6.5–9.0 km)
(Fig. 3). Fast ice flow (>30 m a^{−1}) can be
observed between elevations of 4700 and 4775 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (4.0–5.0 km), where the ice surface velocities during the summer are
approximately 6 m a^{−1} greater than the annual mean
velocity (<40 m a^{−1}). Measurements of winter ice
surface velocities (<10 m a^{−1}) are only available near
the glacier terminus showing a clear interannual variation of the ice
flow speeds.

In August 2009 and 2010, we drilled three 25 m deep shallow
boreholes on LHG12 (Fig. 1b). One borehole was
drilled in the upper ablation area (site 2, approximately
4900 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) and two boreholes were drilled at the automatic weather station (AWS)
locations (sites 1 and 3). The snow and ice temperatures were measured in
the boreholes during the period of 1 October 2010–30 September
2011. The seasonal variations of the snow and ice temperatures in the
shallow boreholes are presented in Figs. 4a, b, and
c. Our measurements show very little fluctuation ($\pm \mathrm{0.4}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$) in the ice temperatures over the depth range
of 20–25 m. Below the 3 m depth, the annual mean
temperature profiles for sites 1 and 2 show a linearly increase in
temperature with depth, while the annual mean temperature profile for
site 3 is convex upward. The mean annual 20 m ice temperatures
(*T*_{20 m}) at sites 1, 2, and 3 are 5.5, 3.0, and
9.5 ^{∘}C higher than the mean annual air temperatures
(*T*_{air}), respectively. Despite its higher elevation, the
near-surface snow and ice temperatures below a depth of 5 m at
site 3 are greater than those in the ablation area (sites 1 and 2),
largely due to the latent heat released as the meltwater entrapped in
the surface snow layers refreezes.

To determine the englacial thermal conditions of LHG12, we drilled an ice core to a depth of 167 m in the upper ablation area of LHG12 (4971 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$; Fig. 1b). In October 2011, the ice temperature was measured to a depth of approximately 110 m using a thermistor string after 20 days of the drilling, as shown in Fig. 4d. The string consists of 50 temperature sensors with a vertical spacing of 0.5 and 10 m at the ice depths of 0–20 and 20–110 m, respectively. The accuracy of the temperature sensor is around $\pm \mathrm{0.05}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$ (Liu et al., 2009). From Fig. 4d we can see that the temperature profile is close to linear with a temperature gradient of around $\mathrm{0.1}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}$ at depths of 9–30 m. Below the depth of 30 m, the ice temperature demonstrates a linear relationship with depth as well but with a smaller temperature gradient of around $\mathrm{0.034}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}$.

Two AWSs were deployed on LHG12, one in the ablation area at 4550 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (site 1, Fig. 1b) and the other in the accumulation area at 5040 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (site 3). During the period of 2010–2013, the mean annual air temperatures (2 m above the ice surface) at sites 1 and 3 were −9.2 and $-\mathrm{12.2}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$, respectively, indicating a lapse rate of $-\mathrm{0.0061}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{1}}$.

The historical knowledge of the surface air temperature and the ELA of LHG12 is a necessity for running the transient model. We reconstruct the past air temperature on LHG12 based on the observations of the Tuole meteorological station (3367 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), which is approximately 175 km southeast to LHG12. We get the precipitation on LHG12 by downscaling the CAPD (China Alpine Region Month Precipitation Dataset) on Qilian Shan. CAPD provides the monthly sum of precipitation during the period of 1960–2013 with a grid spacing of 1 km. We calculate the precipitation on LHG12 from its surrounding 91 grids in CAPD based on the relationship between precipitation and geometric parameters. More details about the reconstruction of both air temperature and precipitation for LHG12 can be found in the Supplement.

3 Thermomechanical ice flow model

Back to toptop
In this study, we use the same 2-D, thermomechanically coupled, first-order flow-band ice flow model as Zhang et al. (2013). Therefore, we only address a very brief review of the model here.

We define *x*, *y*, and *z* as the horizontal along-flow, horizontal
across-flow, and vertical coordinates, respectively. By assuming the
vertical normal stress as hydrostatic and neglecting the bridging
effects (Blatter, 1995; Pattyn, 2002), the equation for momentum
balance is given as

$$\begin{array}{}\text{(1)}& {\displaystyle}{\displaystyle \frac{\partial}{\partial x}}(\mathrm{2}{\mathit{\sigma}}_{xx}^{\prime}+{\mathit{\sigma}}_{yy}^{\prime})+{\displaystyle \frac{\partial {\mathit{\sigma}}_{xy}^{\prime}}{\partial y}}+{\displaystyle \frac{\partial {\mathit{\sigma}}_{xz}^{\prime}}{\partial z}}=\mathit{\rho}g{\displaystyle \frac{\partial s}{\partial x}},\end{array}$$

where ${\mathit{\sigma}}_{ij}^{\prime}$ is the deviatoric stress tensor, *ρ* is the
ice density, *g* is the gravitational acceleration, and *s* is the ice
surface elevation. The parameters used in this study are given in
Table 1.

The constitutive relationship of ice dynamics is described by Glen's flow law (Cuffey and Paterson, 2010):

$$\begin{array}{}\text{(2)}& {\displaystyle}{\displaystyle}{\mathit{\sigma}}_{ij}^{\prime}=\mathrm{2}\mathit{\eta}{\dot{\mathit{\u03f5}}}_{ij},\phantom{\rule{1em}{0ex}}\mathit{\eta}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\mathrm{1}/n}({\dot{\mathit{\u03f5}}}_{e}+{\dot{\mathit{\u03f5}}}_{\mathrm{0}}{)}^{(\mathrm{1}-n)/n},\end{array}$$

where *η* is the ice viscosity, ${\dot{\mathit{\u03f5}}}_{ij}$ is the strain
rate, *n* is the flow law exponent, *A* is the flow rate factor,
${\dot{\mathit{\u03f5}}}_{e}$ is the effective strain rate, and
${\dot{\mathit{\u03f5}}}_{\mathrm{0}}$ is a small number used to avoid singularity. The
flow rate factor is parameterized using the Arrhenius relationship as

$$\begin{array}{}\text{(3)}& {\displaystyle}{\displaystyle}A\left(T\right)={A}_{\mathrm{0}}\mathrm{exp}\left(-{\displaystyle \frac{Q}{RT}}\right),\end{array}$$

where *A*_{0} is the pre-exponential constant, *Q* is the activation
energy for creep, *R* is the universal gas constant, and *T* is the
ice temperature. The effective strain rate ${\dot{\mathit{\u03f5}}}_{e}$ is
related to the velocity gradient by

$$\begin{array}{}\text{(4)}& {\displaystyle}{\displaystyle}{\dot{\mathit{\u03f5}}}_{e}^{\mathrm{2}}\simeq {\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial y}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial y}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial z}}\right)}^{\mathrm{2}},\end{array}$$

where *u* and *v* are the velocity components along the *x* and
*y* direction, respectively. By assuming $\partial v/\partial y=(u/W)(\partial W/\partial x)$, we parameterize the lateral drag,
${\mathit{\sigma}}_{xy}^{\prime}$, as a function of the flow-band half width, *W*,
following Flowers et al. (2011):

$$\begin{array}{}\text{(5)}& {\displaystyle}{\displaystyle}{\mathit{\sigma}}_{xy}^{\prime}=-{\displaystyle \frac{\mathit{\eta}u}{W}}.\end{array}$$

For an easy numerical implementation, we reformulate the momentum balance equation (Eq. 1) as

$$\begin{array}{ll}{\displaystyle \frac{u}{W}}& {\displaystyle}\left\{\mathrm{2}{\displaystyle \frac{\partial \mathit{\eta}}{\partial x}}{\displaystyle \frac{\partial W}{\partial x}}+\mathrm{2}\mathit{\eta}\left[{\displaystyle \frac{{\partial}^{\mathrm{2}}W}{\partial {x}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{W}}{\left({\displaystyle \frac{\partial W}{\partial x}}\right)}^{\mathrm{2}}\right]-{\displaystyle \frac{\mathit{\eta}}{W}}\right\}\\ {\displaystyle}& {\displaystyle}+{\displaystyle \frac{\partial u}{\partial x}}\left(\mathrm{4}{\displaystyle \frac{\partial \mathit{\eta}}{\partial x}}+{\displaystyle \frac{\mathrm{2}\mathit{\eta}}{W}}{\displaystyle \frac{\partial W}{\partial x}}\right)+{\displaystyle \frac{\partial u}{\partial z}}{\displaystyle \frac{\partial \mathit{\eta}}{\partial z}}+\mathrm{4}\mathit{\eta}{\displaystyle \frac{{\partial}^{\mathrm{2}}u}{\partial {x}^{\mathrm{2}}}}+\mathit{\eta}{\displaystyle \frac{{\partial}^{\mathrm{2}}u}{\partial {z}^{\mathrm{2}}}}\\ \text{(6)}& {\displaystyle}& {\displaystyle}=\mathit{\rho}g{\displaystyle \frac{\partial s}{\partial x}},\end{array}$$

where the ice viscosity is expressed as

$$\begin{array}{ll}{\displaystyle}\mathit{\eta}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\mathrm{1}/n}& {\displaystyle}\left[{\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{u}{W}}{\displaystyle \frac{\partial W}{\partial x}}\right)}^{\mathrm{2}}+{\displaystyle \frac{u}{W}}{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial W}{\partial x}}\right.\\ \text{(7)}& {\displaystyle}& {\displaystyle}{\left.+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial z}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{u}{W}}\right)}^{\mathrm{2}}+{{\dot{\mathit{\u03f5}}}_{\mathrm{0}}}^{\mathrm{2}}\right]}^{(\mathrm{1}-n)/\mathrm{2}n}.\end{array}$$

The ice temperature field can be calculated using a 2-D heat transfer equation (Pattyn, 2002):

$$\begin{array}{}\text{(8)}& {\displaystyle}{\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{k}{\mathit{\rho}{c}_{p}}}\left({\displaystyle \frac{{\partial}^{\mathrm{2}}T}{\partial {x}^{\mathrm{2}}}}+{\displaystyle \frac{{\partial}^{\mathrm{2}}T}{\partial {z}^{\mathrm{2}}}}\right)-\left(u{\displaystyle \frac{\partial T}{\partial x}}+w{\displaystyle \frac{\partial T}{\partial z}}\right)+{\displaystyle \frac{\mathrm{4}\mathit{\eta}{{\dot{\mathit{\u03f5}}}_{e}}^{\mathrm{2}}}{\mathit{\rho}{c}_{p}}},\end{array}$$

where *w* is the vertical ice velocity, and *k* and *c*_{p} are the thermal
conductivity and heat capacity of the ice, respectively.

The pressure melting point of the ice, *T*_{pmp}, is described by
the Clausius–Clapeyron relationship

$$\begin{array}{}\text{(9)}& {\displaystyle}{\displaystyle}{T}_{\text{pmp}}={T}_{\mathrm{0}}-\mathit{\beta}(s-z),\end{array}$$

where *T*_{0} is the triple-point temperature of water and *β* is
the Clausius–Clapeyron constant. Following Zhang et al. (2013), we
determined the position of the cold-temperate ice transition surface
(CTS) by considering the following two cases: (1) melting condition,
i.e. cold ice flows downward into the temperate ice zone (TIZ), and (2)
freezing condition, i.e. temperate ice flows upward into the cold ice
zone (Blatter and Greve, 2015; Blatter and Hutter, 1991). For the melting case, the ice
temperature profile at the CTS simply follows a Clausius–Clapeyron
gradient (*β*). However, for the freezing case, the latent heat,
*Q*_{r}, that is released when the water contained in the temperate zone
refreezes is determined as (Funk et al., 1994)

$$\begin{array}{}\text{(10)}& {\displaystyle}{\displaystyle}{Q}_{r}=\mathit{\mu}w{\mathit{\rho}}_{\mathrm{w}}L,\end{array}$$

where *μ* is the fractional water content of the temperate ice,
*ρ*_{w} is the water density, and *L* is the latent heat of
freezing. In this case, following Funk et al. (1994), the ice
temperature gradient at the CTS can be described as

$$\begin{array}{}\text{(11)}& {\displaystyle}{\displaystyle \frac{\partial T}{\partial z}}=-{\displaystyle \frac{{Q}_{r}}{k}}+\mathit{\beta}.\end{array}$$

The free surface evolution follows the kinematic boundary equation

$$\begin{array}{}\text{(12)}& {\displaystyle}{\displaystyle \frac{\partial s}{\partial t}}={b}_{n}+{w}_{\mathrm{s}}-{u}_{\mathrm{s}}{\displaystyle \frac{\partial s}{\partial x}},\end{array}$$

where *s*(*x*,*t*) is the free surface elevation, *u*_{s} and *w*_{s} are the
surface velocity components in *x* and *z*, and *b*_{n} is the surface
mass balance.

For the ice flow model, we assume a stress-free condition for the glacier surface and use the Coulomb friction law to describe the ice–bedrock interface where the ice slips (Schoof, 2005):

$$\begin{array}{}\text{(13)}& {\displaystyle}{\displaystyle}{\mathit{\tau}}_{\mathrm{b}}=\mathrm{\Gamma}{\left({\displaystyle \frac{{u}_{\mathrm{b}}}{{u}_{\mathrm{b}}+{\mathrm{\Gamma}}^{n}{N}^{n}\mathrm{\Lambda}}}\right)}^{\mathrm{1}/n}N,\phantom{\rule{1em}{0ex}}\phantom{\rule{1em}{0ex}}\mathrm{\Lambda}={\displaystyle \frac{{\mathit{\lambda}}_{\text{max}}A}{{m}_{\text{max}}}},\end{array}$$

where *τ*_{b} and *u*_{b} are the basal drag and velocity,
respectively, *N* is the basal effective pressure,
*λ*_{max} is the dominant wavelength of the bed bumps,
*m*_{max} is the maximum slope of the bed bumps, and Γ and
Λ are geometrical parameters (Gagliardini et al., 2007). Here we
take Γ=0.84*m*_{max} following Gagliardini et al. (2007) and
Flowers et al. (2011). The basal effective pressure in the friction law,
*N*, is defined as the difference between the ice overburden pressure
and the basal water pressure (Flowers et al., 2011; Gagliardini et al., 2007):

$$\begin{array}{}\text{(14)}& {\displaystyle}{\displaystyle}N=\mathit{\rho}gH-{P}_{\mathrm{w}}=\mathit{\varphi}\mathit{\rho}gH,\end{array}$$

where *H* and *P*_{w} are the ice thickness and basal water pressure,
respectively, and *ϕ* implies the ratio of basal effective pressure
to the ice overburden pressure. The basal drag is defined as the sum
of all resistive forces (Pattyn, 2002; Van der Veen and Whillans, 1989). It
should be noted that the basal sliding is only permitted when basal
ice temperature reaches the local pressure melting point.

We apply a Dirichlet temperature constraint (*T*_{sbc}) on the
ice surface in the temperature model. In some studies,
*T*_{sbc}=*T*_{air} is used (Zhang et al., 2013),
which, as suggested by recent studies, could result in cold bias in
ice temperature simulations (Meierbachtol et al., 2015). By contrast,
Meierbachtol et al. (2015) recommended using the near-surface
temperature *T*_{dep} at a depth where interannual variations of
air temperatures are damped (15–20 m deep) as a proxy for the
annual mean ice surface temperature. One advantage of using
*T*_{dep} is that the effects of refreezing meltwater and the
thermal insulation of winter snow can be included in the model
(Blatter, 1987; Huang et al., 1982). In fact, the condition
*T*_{dep}=*T*_{air} is acceptable only in dry and cold snow
zones (Cuffey and Paterson, 2010); however,
*T*_{dep}>*T*_{air} is often observed in zones where meltwater
is refreezing in glaciers, such as the LHG12
(Fig. 4). In this study, we set *T*_{sbc} in the
accumulation zone to the glacier near-surface temperature
*T*_{dep}, while *T*_{sbc} in the ablation area is
prescribed by a simple parameterization (Gilbert et al., 2010; Lüthi and Funk, 2001)
as

$$\begin{array}{}\text{(15)}& {\displaystyle}{\displaystyle}{T}_{\text{sbc}}=\left\{\begin{array}{ll}{T}_{\text{dep}},& s>\text{ELA},\\ {T}_{\text{air}}+c,& s\le \text{ELA},\end{array}\right.\end{array}$$

where ELA is the equilibrium line altitude, and *c* is a tuning
parameter implicitly accounting for effects of the surface energy
budget and the lapse rate of air temperature (Gilbert et al., 2010). We
use Eq. (15) as the surface thermal boundary condition of
the reference experiment after comparing another numerical experiment
by setting *T*_{sbc}=*T*_{air} (E-air) (see
Sect. 6.1.2 for details).

At the ice–bedrock interface, we apply the following Neumann-type boundary condition in the temperature model:

$$\begin{array}{}\text{(16)}& {\displaystyle}{\displaystyle \frac{\partial T}{\partial z}}=-{\displaystyle \frac{G}{k}},\end{array}$$

where *G* is the geothermal heat flux. We here use a constant geothermal heat
flux, 40 mW m^{−2}, an in situ measurement from the Dunde ice cap
in the western QS (Huang, 1999), over the entire model domain.

In our model we use a finite difference discretization method and
a terrain-following coordinate transformation. The numerical mesh we
use contains 61 grid points in *x* and 41 layers in *z*. The ice flow
model (Eq. 6) is discretized with a second-order
centered difference scheme while the ice temperature model
(Eq. 8) employs a first-order upstream difference scheme
for the horizontal heat advection term and a node-centered difference
scheme for the vertical heat advection term and the heat diffusion
terms. The velocity and temperature fields are iteratively solved by
a relaxed Picard subspace iteration scheme (De Smedt et al., 2010) in
MATLAB. The free surface evolution (Eq. 12) is solved
using a Crank–Nicolson scheme. More details are given in
Zhang et al. (2013).

4 Surface mass balance model

Back to toptop
In our parameterization of the surface thermal boundary condition,
a transient ELA is important in controlling the extent of accumulation
zone which can be largely warmed by the refreezing of meltwater. We
estimate the annual surface mass balance *b*_{n} of LHG12 during
1960–2013 and determine the position of ELA by *b*_{n}=0.

The daily ablation at elevation *z*, *m*_{a}(*z*), is calculated based on
a degree-day method (Braithwaite and Zhang, 2000):

$$\begin{array}{}\text{(17)}& {\displaystyle}{\displaystyle}{m}_{a}\left(z\right)={f}_{\mathrm{m}}\text{PDD}\left(z\right),\end{array}$$

where *f*_{m} is the degree-day factor, and PDD is the daily positive degree-day sum for glacier surface melt,

$$\begin{array}{}\text{(18)}& {\displaystyle}{\displaystyle}\text{PDD}={T}_{\text{mean}}-{T}_{\mathrm{m}},\end{array}$$

where *T*_{mean} is the daily mean air temperature, and *T*_{m} is
the threshold temperature when melt occurs. Note that surface melt may
happen even if ${T}_{\text{mean}}<\mathrm{0}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$ due to the
positive air temperature during the day, indicating a negative *T*_{m}
in such cases.

As the CAPD dataset only provides monthly precipitation sums, we
estimate the daily precipitation on LHG12 by assuming a uniformly
distributed precipitation over a month. The daily accumulation at
elevation *z*, *m*_{c}(*z*), is calculated as follows:

$$\begin{array}{}\text{(19)}& {\displaystyle}{\displaystyle}{m}_{c}\left(z\right)=\left\{\begin{array}{ll}{f}_{\mathrm{P}}\cdot {P}_{\text{total}},& {T}_{\text{mean}}<{T}_{\text{crit1}},\\ {f}_{\mathrm{P}}\cdot \frac{{T}_{\text{crit2}}-{T}_{\text{mean}}}{{T}_{\text{crit2}}-{T}_{\text{crit1}}}\cdot {P}_{\text{total}},& {T}_{\text{crit1}}\le {T}_{\text{mean}}\le {T}_{\text{crit2}},\\ \mathrm{0},& {T}_{\text{mean}}>{T}_{\text{crit2}},\end{array}\right.\end{array}$$

where *P*_{total} is the daily precipitation, *T*_{crit1}
and *T*_{crit2} are the threshold temperatures for the snow and
rain transition, and *f*_{P} is a tuning parameter for the precipitation
to account for the uncertainties of the gridded CAPD data and the
downscaling method. The model is well calibrated by investigating the
sensitivities of model parameters and by comparing the simulated
results with the observed mass balance during 2010–12 (see
Figs. S5–8 in the Supplement).

5 Simulation strategies

Back to toptop
In order to remove the uncertainties remained in the model initial conditions (including initial topography and model parameters) (Seroussi et al., 2011; Zwinger and Moore, 2009), we allow the free surface to relax for a period of 3 years under a constant present-day surface mass balance and zero basal sliding (see details in the Supplement). The time step for the relaxation experiment is set to 0.1 year. We apply surface relaxation before running all of the diagnostic and transient simulations in this study.

First, we explore the model sensitivities to geometrical bed
parameters (*λ*_{max} and *m*_{max}), ratio of basal
effective pressure (*ϕ*), water content (*μ*), geothermal heat
flux (*G*), and valley shape index (*b*) using diagnostic simulations
with relaxed present-day geometry of LHG12. Then, we calibrate the
parameters in surface thermal boundary condition (ELA_{0},
*T*_{dep}, and *c*), where ELA_{0} is the initial ELA, by first
running a diagnostic simulation with a given set of parameters
(ELA_{0}, *T*_{dep}, and *c*) and then performing a transient
simulation using the diagnostic run as an initial condition with
time-dependent ELA and surface air temperature (see
Sect. 5.3). This numerical process of finding the
optimal parameter set of ELA_{0}, *T*_{dep}, and *c* is repeated
until we find a good agreement between the modelled and measured deep
borehole temperatures in 2011. Finally, we perform experiments to
investigate the sensitivities of the thermomechanical model to
different surface thermal boundary conditions (E-ref-D and E-air). We
also perform four experiments (E-advZ, E-advX, E-strain, and E-slip)
to explore the effects of heat advection, strain heating, and basal
sliding on the thermal distribution and flow dynamic behaviours of
LHG12.

To investigate the impacts of historical climate conditions on the
thermal regime of LHG12, time-dependent simulations are performed from
1961 to 2011. The diagnostic run with the calibrated surface thermal
parameters is used as the initial condition of the transient
simulations. In our simulations, we assume that the surface
temperature (*T*_{dep}) in the accumulation zone is constant in
time and space. Due to a lack of in situ observations (e.g. firn
thickness, firn densities) and coupled heat–water transfer model
(e.g. Gilbert et al., 2012; Wilson and Flowers, 2013), we do not simulate the complex
processes of heat exchange in the accumulation zone. To understand how
the thermal status varies over time in the deep borehole, we design
three transient simulations (E-cold, E-warm, E-ref-T) by setting
*T*_{dep} to −5, −1, and $-\mathrm{1.8}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$ (as
calibrated in our diagnostic simulations; see
Sect. 6.1.1), respectively. We also design two
experiments (E-high and E-low) to explore the impacts of the initial
ELA (ELA_{0}) on the thermal regime of LHG12.

In the transient model, the mean annual air temperature
(*T*_{air}) and ELA are allowed to vary in time. However, we keep
the glacier geometry fixed in the transient simulations for two main
reasons: (1) the tributaries of LHG12, which our flow-band model
neglects, may have non-negligible inflow ice fluxes that impact the
mass continuity equation; and (2) the mean surface elevation change above
4600 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (the confluence area) over 1957–2014 is close
to negligible (approximately −10.4 m, around 4.4 %
of the mean ice thickness) (see Fig. S10 in the Supplement).

6 Simulation results and discussions

Back to toptop
As shown in Figs. 5 and 6, we conduct
a series of sensitivity experiments to investigate the relative
importance of different model parameters (*λ*_{max},
*m*_{max}, *ϕ*, *μ*, *G*, *b*) on ice flow speeds and
TIZ sizes by varying the value of one parameter
while holding the others fixed. We set ELA to 4980 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$,
as observed in the 2011 Landsat image of LHG12. *T*_{dep} is set
to $-\mathrm{2.7}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$, which is the 20 m deep ice
temperature at site 3. *c* is set to $-\mathrm{4.3}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$, which
is the mean differences between the 20 m deep ice temperatures
and mean annual air temperatures at sites 1 and 2.

The friction law parameters, *λ*_{max} and *m*_{max},
which describe the geometries of bedrock obstacles (Flowers et al., 2011; Gagliardini et al., 2007), have non-negligible impacts on
the model results. As shown in Fig. 5a and b, the
modelled velocities and TIZ sizes increase as *λ*_{max}
increases and *m*_{max} decreases, similar to the results
observed by Flowers et al. (2011) and Zhang et al. (2013). A large
increase in the modelled basal sliding velocity occurs when
*m*_{max} is lower than 0.2. The ratio, *ϕ*, is an insensitive
parameter in our model when it is larger than 0.3
(Fig. 5c). Although the water content, *μ*, in the
ice does not directly impact the ice velocity simulations (the flow
rate factor *A* is assumed independent of the water content in ice),
it can affect the temperature field, and consequently influence *A*
and the ice velocities
(Fig. 5d). Figure 6d shows that
increasing the water content may result in larger TIZ size. In
addition, we test the sensitivity of the model to different geothermal
heat flux values. A larger geothermal heat flux can result in a larger
TIZ but has a limited impact on modelled ice velocity
(Figs. 5e and 6e). As shown by
Zhang et al. (2013), our model results are mainly controlled by the
shape of the glacial valley, specifically the *b* index
(Sect. 2.1). A large value of *b* indicates a flat
glacial valley and suggests that a small lateral drag is exerted on
the ice flow (Figs. 5f and 6f).

Based on the sensitivity experiments described above, we adopt
a parameter set of *λ*_{max}=4 m,
*m*_{max}=0.3, *ϕ*=1 (no basal water pressure),
*μ*=3 %, *G*=40 mW m^{−2}, and *b*=1.2 as
a diagnostic reference in our modelling experiment (E-ref-D). To
calibrate the surface thermal parameters, we perform a series of model
runs by varying ELA_{0}, *T*_{dep}, and *c* in the range of
4940–5010 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, −3.3–−1.5, and
1–6 ^{∘}C, respectively. The performance of different
parameter combinations is evaluated by comparing the root mean square
(RMS) of the differences between the measured and modelled temperatures
below the 40 m depth of the deep borehole
(Fig. 7). We find the RMSs are sufficiently
small ($<\mathrm{0.11}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$) when
ELA_{0}=4940 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, ${T}_{\text{dep}}=-\mathrm{1.8}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$, and $c=\mathrm{3}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$
(Fig. 7a–f and h). The sensitivity experiments
of ELA_{0} and their impacts on the thermal regime of LHG12 are
discussed in Sect. 6.2.3.

To investigate the impacts of different surface thermal boundary
conditions on the thermomechanical fields of LHG12, we perform an
experiment (E-air) in which we set *T*_{sbc}=*T*_{air} and
compare its results with those of
E-ref-D. Figure 8a shows that the ice temperatures along
the CL are highly sensitive to *T*_{sbc}. From E-air, it is
observed that LHG12 becomes fully cold (Fig. 8b), with
an average field temperature 6.3 ^{∘}C lower than that
of E-ref-D (Fig. 8a), which decreases the ice surface
velocity by approximately 11 m a^{−1}
(Fig. 8c).

As noted above, the dynamics of LHG12 can be strongly influenced by
the choices of different surface thermal boundary conditions. For
LHG12, most accumulation and ablation events overlap in summer
(Sun et al., 2012). The meltwater entrapped in snow and moulins can
refreeze at a depth where the temperature is below the melting
point. Refreezing of meltwater releases large amounts of latent heat,
which can significantly raise the near-surface snow and ice temperatures
and result in the warm bias of *T*_{20m} compared to
*T*_{air} (Fig. 4). Therefore, compared with
E-air, the experiment E-ref-D better incorporates the effects of
meltwater refreezing in the accumulation basin into the prescribed
surface thermal boundary constraints, resulting in more accurate
simulations of ice temperature and flow fields.

To assess the relative contributions of heat advection and strain
heating to the thermomechanical field of LHG12, we conduct three
experiments (E-advZ, E-advX, and E-strain), in which the vertical
advection, horizontal along-flow advection, and strain heating were
neglected, respectively. In addition, to investigate the effects of
basal sliding predicted by the Coulomb friction law on the thermal
state and flow dynamics of LHG12, we perform an experiment (E-slip)
with *u*_{b}=0.

Figure 8 compares the ice velocity and temperature
results of E-advZ, E-advX, E-strain, and E-slip with those of
E-ref-D. If the vertical advection is neglected (E-advZ), cold ice at
the glacier surface cannot be transported downwards into the interior
of LHG12, and therefore LHG12 becomes warmer
(Fig. 8a and b) and flows faster relative to other
experiments (Fig. 8c). Figure 8a shows
a transition of the mean column ice temperature at 1.8 km (the
horizontal position of ELA) in E-advX. This transition arises from the
discontinuous surface thermal boundary conditions across ELA. Compared
with the reference experiment, E-advX predicts lower field
temperatures (by 2.3 ^{∘}C on average;
Fig. 8a) and much smaller surface ice velocities (by
7.2 m a^{−1} on average; Fig. 8c). Because the
accumulation basin of LHG12 is relatively warm, E-advX, which neglects
the horizontal transport of ice from upstream to downstream, predicts
much colder conditions for LHG12 and the modelled temperate ice only
appears in the accumulation zone (Fig. 8d). We observe
that strain heating contributes greatly to the thermal configuration
of LHG12. If we leave away the strain heating, the modelled mean ice
temperature field is lower than that of E-ref-D by
0.4 ^{∘}C on average (Fig. 8a). The TIZ
becomes much thinner compared to E-ref-D, indicating the importance of
strain heating in the formation of basal temperate ice
(Fig. 8d). Previous studies have suggested that basal
sliding can significantly influence the thermal structures and
velocity fields of glaciers
(Wilson et al., 2013; Zhang et al., 2015). However, in this study
the neglect of basal sliding (E-slip) results in a temperature field
very similar to that of E-ref-D. We attribute this difference to the
relatively small modelled basal sliding values for LHG12.

In the reference experiment (E-ref-T), we simulate the distributions
of horizontal ice velocities and temperatures in a transient manner
(Fig. 9a and c). Next, the model results are
compared to the measured ice surface velocities and the ice
temperature profile in the deep borehole
(Figs. 9b and d). Generally, the modelled ice
surface velocities are in good agreement with in situ observations
from the glacier head to 4.8 km along the CL
(Fig. 9b). However, from 4.8 km to the glacier
terminus, our model generally underestimates the ice surface
velocities, as shown in all simulations in Fig. 5. There
are three possible reasons for this underestimation. First, the model
neglects the convergent ice fluxes from the west branch. Second, an
enhanced basal sliding zone may exist at the confluence area, which is
not captured by the model. Third, the transient model uses a fixed
topography, which may fail to capture the time-dependent glacier
changes that can largely influence the ice flow dynamics. Here we
verify the first and second hypotheses by conducting two other
experiments, i.e. E-W and E-WS. In E-W the glacier widths are
increased by 450 m at 5.8–7.3 km as a proxy of including the
impact of the convergent flow from the west branch
(Fig. 10). In E-WS, except for the same glacier width
increasing as in E-W, we also increase *λ*_{max} by
100 % and decrease *m*_{max} by 33 % for
accelerating the basal sliding at 5.8–7.3 km
(Fig. 10). We can clearly find that while both factors have
a non-negligible contribution to the model results, the basal sliding
may play a bit more important role in the confluence area. This
indicates a need of considering glacier flow branches and spatially
variable sliding law parameters in real glacier modelling studies.

The model predicts a TIZ overlain by cold ice over a horizontal
distance of 0.6–6.5 km (Fig. 9c). In addition,
we further compare our model results with the in situ 110 m
deep ice temperature measurements
(Fig. 9d). Modelled and measured borehole
temperature profiles show a close match within a RMS
difference of 0.2 ^{∘}C below the 20 m
depth. Because in situ ice temperature data below the 110 m
depth have not been obtained, we are unable to compare the modelled and
measured ice temperatures at the ice–bedrock interface.

The transient temperature profiles in our reference experiment
(E-ref-T) show a cooling trend in the upper surface (0–50 m
deep) particularly during the period 1991–2011
(Fig. 11a). This surface cooling trend is also
suggested by the experiments E-cold and E-warm in which the thermal
parameter *T*_{dep} is set to −5 and $-\mathrm{1}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$,
respectively. However, E-cold (E-warm) underestimates (overestimates)
the ice temperature in the deep borehole compared with the
observation. From temperature profile modelled by the reference
diagnostic experiment (E-ref-D; Fig. 11a), we can
see that the shape of the upper temperature profile (0–30 m)
is still difficult to simulate.

In the surface thermal boundary condition, the parameter
*T*_{dep} simply represent a mean thermal status of the
near-surface region in the accumulation zone. The reference transient
experiment, which adopts the tuned value of *T*_{dep}
($-\mathrm{1.8}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$), shows a good agreement between the
modelled and observed temperature profile, suggesting a possibly
supportive evidence that the calibrated *T*_{dep} may indeed
contain part of the historical climate information of LHG12. In
Fig. 11b, we can see that both the summer (June,
July and August) air temperature and ELA appear a slight (large)
increase during 1971–1991 (1991–2011), which can explain the small
(large) decrease of ice temperatures for all model experiments over
the same time period, indicating that LHG12 (or similar type glaciers)
may not accordingly become cold under a cold climate scenario, when
a relatively large accumulation basin grows and more refreezing latent
heat from meltwater is released.

For LHG12, ELA determines the size of accumulation zone and the amount
of latent heat released by meltwater refreezing. In the above
transient simulations (see Sects. 6.2.1 and
6.2.2), the initial conditions are generated by
diagnostic runs in which the ELA_{0} is set to 4940 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$
To investigate the impacts of initial ELA on the thermal regime of
LHG12, we perform two additional experiments with different initial
ELAs, i.e. E-high (ELA_{0}=5000 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) and E-low
(ELA_{0}=4800 $\mathrm{m}\phantom{\rule{0.125em}{0ex}}\mathrm{a}.\mathrm{s}.\mathrm{l}.$).

The initial differences of the column mean temperature between the
sensitivity experiments (E-high and E-low) and the reference
experiment (E-ref-T) are significant near the downstream of the
ELA_{0}, as can be seen in Fig. 12. After 50-year
runs, however, they are largely reduced because we set the same
surface thermal boundary conditions for E-high/E-low and E-ref-T
during the transient simulations. However, we can still clearly observe two
phase shifts of temperature difference in space between the model
initialization and the year 2011. By heat advection the temperature
field of the downstream and the deep part of LHG12 can still be
remarkably impacted by the thermal status of the upstream after dozens
of years, which suggests that the reconstruction of past climate
changes on the glacier is crucial to better estimate the change of
englacial thermal condition (Gilbert et al., 2014a).

Although our 2-D, first-order flow-band model can account for part of
the three-dimensional nature of LHG12 by parameterizing the lateral
drag with glacier width variations, it cannot fully describe the ice
flow along the *y* direction and is not able to account for the
confluence of glacier tributaries. The shape of the LHG12 glacier
valley is described using a constant value for index *b* (1.2;
approximately V-type cross sections), which was determined from
several traverse GPR profiles (Fig. 1). However, for
real glaciers, the cross-sectional geometry profiles are generally
complex, resulting in an inevitable bias when we idealize the glacier
cross-sectional profiles by using power law functions across the
entire LHG12. Although the regularized Coulomb friction law provides
a physical relationship between the basal drag and sliding velocities,
several parameters (e.g. *λ*_{max}, *m*_{max}) still
must be prescribed based on surface velocity observations. Another
uncertainty could be from the spatially uniform geothermal heat flux
that we assume in the model, as it may have a great spatial variation
due to mountain topography (Lüthi and Funk, 2001). In addition, we can
also improve our model ability by linking the water content in the
temperate ice layer to a physical thermo-hydrological process in the
future.

Due to the limitations of in situ shallow borehole ice temperature measurements, the surface thermal boundary condition in our temperature model is determined using a simple parameterization based on observations at three elevations (Fig. 1b). In addition, the parameterized surface thermal boundary condition only provides a rough estimate of the overall contributions of the heat from refreezing meltwater and ice flow advection. At this stage, we cannot simulate the actual physical process involved in the transport of near-surface heat from refreezing, which has been suggested by Gilbert et al. (2012) and Wilson and Flowers (2013). In our transient simulations, the glacier geometry and surface thermal parameters are assumed constant in time, which may lead to some unphysical model outputs. Due to a lack of long-term historical climate and geometry inputs, our transient simulations are only performed during 1961–2011 initialized with according diagnostic runs. However, the transient simulations in this study are mainly for seeking some possible reasons for the formation of the current temperature profile of the deep borehole. Therefore, we do not expect accurate model results from the transient experiments.

7 Conclusions

Back to toptop
For the first time, we investigate the thermomechanical features of a typical valley glacier, Laohugou Glacier No. 12, on Qilian Shan, which is an important fresh water source for the arid regions in western China. We assess the thermomechanical features of LHG12 using a two-dimensional, thermomechanically coupled, first-order flow-band model based on available in situ measurements of glacier geometries, borehole ice temperatures, and surface meteorological and velocity observations.

Similar to other alpine land-terminating glaciers, the mean annual
horizontal ice flow speeds of LHG12 are relatively low (less than
40 m a^{−1}). However, we observed large interannual
variations in the ice surface velocity during the summer and winter
seasons. Due to the release of heat from refreezing meltwater, the
observed ice temperatures for the shallow ice borehole in the
accumulation basin (site 3; Fig. 1) are higher than
those at sites 1 and 2 at lower elevations, indicating the existence
of meltwater refreezing. Thus, we parameterize the surface thermal
boundary condition by accounting for the 20 m deep temperature
instead of only the surface air temperatures. We observed that LHG12
has a polythermal structure with a TIZ that is overlain
by cold ice near the glacier base throughout a large region of the
ablation area. Time-dependent simulations reveal that the englacial
temperature becomes colder in recent two decades as a consequence of
the shrink of accumulation area and rising surface air temperature.

Horizontal heat advection is important on LHG12 for bringing the
relatively warm ice in the accumulation basin (due to the heat from
refreezing meltwater) to the downstream ablation zone. In addition,
vertical heat advection is important for transporting the near-surface
cold ice downwards into the glacier interior, which “cools down” the
ice. Furthermore, we argue that the strain heating of LHG12 also plays
an important role in controlling the englacial thermal status, as
suggested by Zhang et al. (2015).
However, we also observed that
simulated basal sliding (very small; <4 m a^{−1})
contributes little to the thermomechanical configuration of LHG12.

The mean annual surface air temperature could serve as a good approximation for the temperatures of shallow ice, where seasonal climate variations are damped at cold and dry locations (Cuffey and Paterson, 2010). However, for LHG12, using the mean annual surface air temperature as the thermal boundary condition at the ice surface would predict an entirely cold glacier with very small ice flow speeds. For LHG12, a decline of ELA under a cold climate may assist an increase of the amount of refreezing latent heat in the accumulation basin and, therefore, possibly raise the englacial temperature. Because warming is occurring on alpine glaciers on, for example, the Himalayas and Qilian Shan, further studies of supraglacial and near-surface heat transport are very important because they will affect the surface thermal conditions and, eventually, the dynamical behaviours of the glacier.

Data availability

Back to toptop
Data availability.

All the input data, results, and codes are available upon request by email to Yuzhe Wang (wangyuzhe@lzb.ac.cn).

Competing interests

Back to toptop
Competing interests.

The authors declare that they have no conflict of interest.

Supplement

Back to toptop
Supplement.

The supplement related to this article is available online at: https://doi.org/10.5194/tc-12-851-2018-supplement.

Acknowledgements

Back to toptop
Acknowledgements.

This work is supported by the National Basic Research Program (973) of
China (2013CBA01801, 2013CBA01804) and the Technology Services
Network Program (STS-HHS Program) of Cold and Arid Regions
Environmental and Engineering Research Institute, Chinese Academy of
Sciences (HHS-TSS-STS-1501). Tong Zhang is also supported by the
National Natural Science Foundation of China (41601070). We are
grateful to numerous people for their hard fieldwork and for the
support from the Qilian Shan Station of Glaciology and Ecologic
Environment, Chinese Academy of Sciences. We thank Chen Rensheng
and Liu Junfeng for providing the CAPD precipitation dataset. We thank
Martin Lüthi, Andy Aschwanden, and an anonymous reviewer for their
thorough and constructive reviews which significantly improved the
manuscript, as well as Oliver Gagliardini for his editorial
work.

Edited by: Olivier
Gagliardini

Reviewed by: Martin Lüthi, Andy
Aschwanden, and one anonymous referee

References

Back to toptop
Blatter, H.: On the thermal regime of an Arctic valley glacier: a study of White Glacier, Axel Heiberg Island, NWT, Canada, J. Glaciol., 33, 200–211, https://doi.org/10.3189/S0022143000008704, 1987. a, b

Blatter, H.: Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients, J. Glaciol., 41, 333–344, 1995. a

Blatter, H. and Greve, R.: Comparison and verification of enthalpy schemes for polythermal glaciers and ice sheets with a one-dimensional model, Polar Sci., 9, 196–207, 2015. a

Blatter, H. and Hutter, K.: Polythermal conditions in Arctic glaciers, J. Glaciol., 37, 261–269, 1991. a, b

Blatter, H. and Kappenberger, G.: Mass balance and thermal regime of Laika ice cap, Coburg Island, NWT, Canada, J. Glaciol., 34, 102–110, https://doi.org/10.3189/S0022143000009126, 1988. a

Braithwaite, R. J. and Zhang, Y.: Sensitivity of mass balance of five Swiss glaciers to temperature changes assessed by tuning a degree-day model, J. Glaciol., 46, 7–14, https://doi.org/10.3189/172756500781833511, 2000. a

Cuffey, K. and Paterson, W.: The Physics of Glaciers, 4th edn., Elsevier, Amsterdam, 2010. a, b, c

De Smedt, B., Pattyn, F., and De Groen, P.: Using the unstable manifold correction in a Picard iteration to solve the velocity field in higher-order ice-flow models, J. Glaciol., 56, 257–261, 2010. a

Du, W., Qin, X., Liu, Y., and Wang, X.: Variation of the Laohugou Glacier No. 12 in the Qilian Mountains, Journal of Glaciology and Geocryology, 30, 373–379, 2008 (in Chinese with English summary). a

Flowers, G. E., Roux, N., Pimentel, S., and Schoof, C. G.: Present dynamics and future prognosis of a slowly surging glacier, The Cryosphere, 5, 299–313, https://doi.org/10.5194/tc-5-299-2011, 2011. a, b, c, d, e

Funk, M., Echelmeyer, K. A., and Iken, A.: Mechanisms of fast flow in Jakobshavn Isbræ, West Greenland: Part II. Modeling of englacial temperatures, J. Glaciol., 40, 569–585, 1994. a, b

Gagliardini, O., Cohen, D., Råback, P., and Zwinger, T.: Finite-element modeling of subglacial cavities and related friction law, J. Geophys. Res., 112, F02027, https://doi.org/10.1029/2006JF000576, 2007. a, b, c, d

Gilbert, A., Wagnon, P., Vincent, C., Ginot, P., and Funk, M.: Atmospheric warming at a high-elevation tropical site revealed by englacial temperatures at Illimani, Bolivia (6340 m above sea level, 16^{∘} S, 67^{∘} W), J. Geophys. Res., 115, D10109, https://doi.org/10.1029/2009JD012961, 2010. a, b, c

Gilbert, A., Vincent, C., Wagnon, P., Thibert, E., and Rabatel, A.: The influence of snow cover thickness on the thermal regime of Tête Rousse Glacier (Mont Blanc range, 3200 m a.s.l.): Consequences for outburst flood hazards and glacier response to climate change, J. Geophys. Res., 117, F04018, https://doi.org/10.1029/2011JF002258, 2012. a, b, c

Gilbert, A., Gagliardini, O., Vincent, C., and Wagnon, P.: A 3-D thermal regime model suitable for cold accumulation zones of polythermal mountain glaciers, J. Geophys. Res.-Earth, 119, 1876–1893, 2014a. a, b

Gilbert, A., Vincent, C., Six, D., Wagnon, P., Piard, L., and Ginot, P.: Modeling near-surface firn temperature in a cold accumulation zone (Col du Dôme, French Alps): from a physical to a semi-parameterized approach, The Cryosphere, 8, 689–703, https://doi.org/10.5194/tc-8-689-2014, 2014b. a

Gilbert, A., Vincent, C., Gagliardini, O., Krug, J., and Berthier, E.: Assessment of thermal change in cold avalanching glaciers in relation to climate warming, Geophys. Res. Lett., 42, 6382–6390, 2015. a

Guo, W., Xu, J., Liu, S., Shangguan, D., Wu, L., Yao, X., Zhao, J., Liu, Q., Jiang, Z., Li, P., Wei, J., Bao, W., Yu, P., Ding, L., Li, G., Ge, C., and Wang, Y.: The second glacier inventory dataset of China (Version 1.0), Cold and Arid Regions Science Data Center at Lanzhou, https://doi.org/10.3972/glacier.001.2013.db, 2014. a

Guo, W., Liu, S., Xu, J., Wu, L., Shangguan, D., Yao, X., Wei, J., Bao, W., Yu, P., Liu, Q., and Jiang, Z.: The second Chinese glacier inventory: data, methods and results, J. Glaciol., 61, 357–372, 2015. a

Huang, M.: On the temperature distribution of glaciers in China, J. Glaciol., 36, 210–216, 1990. a

Huang, M.: Forty year's study of glacier temperature in China, Journal of Glaciology and Geocryology, 21, 193–199, 1999 (in Chinese with English summary). a

Huang, M., Wang, Z., and Ren, J.: On the temperature regime of continental-type glaciers in China, J. Glaciol., 28, 117–128, https://doi.org/10.3189/S0022143000011837, 1982. a, b

Huss, M., Jouvet, G., Farinotti, D., and Bauder, A.: Future high-mountain hydrology: a new parameterization of glacier retreat, Hydrol. Earth Syst. Sci., 14, 815–829, https://doi.org/10.5194/hess-14-815-2010, 2010. a

Irvine-Fynn, T. D. L., Hodson, A. J., Moorman, B. J., Vatne, G., and Hubbard, A. L.: Polythermal glacier hydrology: a review, Rev. Geophys., 49, RG4002, https://doi.org/10.1029/2010RG000350, 2011. a

Jansson, P., Hock, R., and Schneider, T.: The concept of glacier storage: a review, J. Hydrol., 282, 116–129, 2003. a

Liu, Y., Hou, S., Wang, Y., and Song, L.: Distribution of borehole temperature at four high-altitude alpine glaciers in Central Asia, J. Mt. Sci., 6, 221–227, 2009. a

Liu, Y., Qin, X., Du, W., Sun, W., and Hou, D.: The movement features analysis of Laohugou Glacier No. 12 in Qilian Mountains, Sciences in Cold and Arid Regions, 3, 119–123, 2011. a, b

Lüthi, M. P. and Funk, M.: Modelling heat flow in a cold, high-altitude glacier: interpretation of measurements from Colle Gnifetti, Swiss Alps, J. Glaciol., 47, 314–324, 2001. a, b

Meierbachtol, T. W., Harper, J. T., Johnson, J. V., Humphrey, N. F., and Brinkerhoff, D. J.: Thermal boundary conditions on western Greenland: observational constraints and impacts on the modeled thermomechanical state, J. Geophys. Res.-Earth, 120, 623–636, 2015. a, b, c, d

Müller, F.: On the thermal regime of a High-Arctic valley glacier, J. Glaciol., 16, 119–133, https://doi.org/10.3189/S0022143000031476, 1976. a

Neckel, N., Kropáček, J., Bolch, T., and Hochschild, V.: Glacier mass changes on the Tibetan Plateau 2003–2009 derived from ICESat laser altimetry measurements, Environ. Res. Lett., 9, 014009, 2014. a

Pattyn, F.: Transient glacier response with a higher-order numerical ice-flow model, J. Glaciol., 48, 467–477, 2002. a, b, c

Schäfer, M., Gillet-Chaulet, F., Gladstone, R., Pettersson, R., A. Pohjola, V., Strozzi, T., and Zwinger, T.: Assessment of heat sources on the control of fast flow of Vestfonna ice cap, Svalbard, The Cryosphere, 8, 1951–1973, https://doi.org/10.5194/tc-8-1951-2014, 2014. a

Schoof, C.: The effect of cavitation on glacier sliding, P. R. Soc. A, 461, 609–627, 2005. a

Seroussi, H., Morlighem, M., Rignot, E., Larour, E., Aubry, D., Ben Dhia, H., and Kristensen, S. S.: Ice flux divergence anomalies on 79north Glacier, Greenland, Geophys. Res. Lett., 38, https://doi.org/10.1029/2011GL047338, 2011. a

Shi, Y. and Liu, S.: Estimation on the response of glaciers in China to the global warming in the 21st century, Chinese Sci. Bull., 45, 668–672, 2000. a

Sun, W., Qin, X., Ren, J., Yang, X., Zhang, T., Liu, Y., Cui, X., and Du, W.: The surface energy budget in the accumulation zone of the Laohugou Glacier No. 12 in the western Qilian Mountains, China, in summer 2009, Arct. Antarct. Alp. Res., 44, 296–305, 2012. a

Svensson, H.: Is the cross-section of a glacial valley a parabola?, J. Glaciol., 3, 362–363, 1959. a

Tian, H., Yang, T., and Liu, Q.: Climate change and glacier area shrinkage in the Qilian Mountains, China, from 1956 to 2010, Ann. Glaciol., 55, 187–197, 2014. a

Van der Veen, C. and Whillans, I.: Force budget: I. Theory and numerical methods, J. Glaciol., 35, 53–60, 1989. a

Vincent, C., Le Meur, E., Six, D., Possenti, P., Lefebvre, E., and Funk, M.: Climate warming revealed by englacial temperatures at Col du Dôme (4250 m, Mont Blanc area), Geophys. Res. Lett., 34, L16502, https://doi.org/10.1029/2007GL029933, 2007. a

Wang, Y., Ren, J., Qin, X., Liu, Y., Zhang, T., Chen, J., Li, Y., and Qin, D.: Ice depth and glacier-bed characteristics of Laohugou Glacier No. 12, Qilian Mountains, revealed by ground-penetrating radar, Journal of Glaciology and Geocryology, 1, 28–35, 2016 (in Chinese with English summary). a, b

Wilson, N. J. and Flowers, G. E.: Environmental controls on the thermal structure of alpine glaciers, The Cryosphere, 7, 167–182, https://doi.org/10.5194/tc-7-167-2013, 2013. a, b, c

Wilson, N. J., Flowers, G. E., and Mingo, L.: Comparison of thermal structure and evolution between neighboring subarctic glaciers, J. Geophys. Res.-Earth, 118, 1443–1459, 2013. a, b

Yao, T., Thompson, L., Yang, W., Yu, W., Gao, Y., Guo, X., Yang, X., Duan, K., Zhao, H., and Xu, B.: Different glacier status with atmospheric circulations in Tibetan Plateau and surroundings, Nat. Clim. Change, 2, 663–667, 2012. a

Zhang, T., Xiao, C., Colgan, W., Qin, X., Du, W., Sun, W., Liu, Y., and Ding, M.: Observed and modelled ice temperature and velocity along the main flowline of East Rongbuk Glacier, Qomolangma (Mount Everest), Himalaya, J. Glaciol., 59, 438–448, 2013. a, b, c, d, e, f, g

Zhang, T., Ju, L., Leng, W., Price, S., and Gunzburger, M.: Thermomechanically coupled modelling for land-terminating glaciers: a comparison of two-dimensional, first-order and three-dimensional, full-Stokes approaches, J. Glaciol., 61, 702–712, 2015. a, b

Zwinger, T. and Moore, J. C.: Diagnostic and prognostic simulations with a full Stokes model accounting for superimposed ice of Midtre Lovénbreen, Svalbard, The Cryosphere, 3, 217–229, https://doi.org/10.5194/tc-3-217-2009, 2009. a