A mathematical model of polythermal glaciers and ice sheets

Abstract

A mathematical model of the flow and temperature distribution of polythermal glaciers or ice sheets is deduced. Cold ice is treated as a non-linear viscous heat conducting fluid, while temperate ice is regarded as a binary mixture of ice and water. The simplest mixture concept with two balance laws of mass but only one balance law of momentum and energy is proposed. The field equations for the ice and water content and the boundary conditions which must hold at the free surface, at the ice-water interface, at the cold-temperature transition surface and at the rock-bed are deduced. In particular it is shown that an earlier formulation of polythermal ice due to Fowler and Larson (1978) is inconsistent. No boundary value problems are solved as the emphasis is on the physical motivation and justification of the principles.     

Conditions for a steady ice sheet–ice shelf junction

This paper investigates the conditions under which a marine ice sheet may adopt a steady profile. The ice is treated as a linear viscous fluid caused to flow from a rigid base to and over water, treated as a denser but inviscid fluid. The solutions in the region around the point of flotation, or ‘transition’ zone, are calculated numerically. In-flow and out-flow conditions appropriate to ice sheet and ice shelf flow are applied at the ends of the transition zone and the rigid base is specified; the flow and steady free surfaces are determined as part of the solutions. The basal stress upstream, and the basal deflection downstream, of the flotation point are examined to determine which of these steady solutions satisfy ‘contact’ conditions that would prevent (i) the steady downstream basal deflection contacting the downstream base, and (ii) the upstream ice commencing to float in the event it was melted at the base. In the case that the upstream bed is allowed to slide, we find only one mass flux that satisfies the contact conditions. When no sliding is allowed at the bed, however, we find a range of mass fluxes satisfy the contact conditions. The effect of ‘backpressure’ on the solutions is investigated, and is found to have no affect on the qualitative behaviour of the junctions. To the extent that the numerical, linearly viscous treatment may be applied to the case of ice flowing out over the ocean, we conclude that when sliding is present, Weertman's ‘instability’ hypothesis holds.     

On similarity solution of an unsteady laminar boundary layer along a flat plate

Summary A unified analysis has been made to obtain all possible similarity solutions of the steady and unsteady, forced flow, inside a boundary layer along a flat plate. Though previously, attempts were made to obtain similarity solutions of a steady boundary layer flow neglecting viscous dissipation term in the energy conservation equation but the treatments were not complete. Here we have taken account of the viscous dissipation term. In the steady case it has been shown that for a similarity solution of both velocity and temperature, there should be a relation between the undisturbed flow outside the boundary layer and the temperature of the plate. It has been shown that the similarity solution exists in the unsteady case if we neglet the viscous dissipation term in the energy equation.     

Convection and gravity waves in two layer models. I. Overstable modes driven in conducting boundary layers

Abstract

Stability analysis is formulated for a two-layer fluid model in which the upper and lower layers are convectively stable and unstable, respectively. With discontinuities in viscosity and conductivity at the interface, the exchange of stability does not generally hold and overstability is possible. A detailed analytical treatment is presented for the case of small viscosity and conductivity in which viscous and conducting boundary layers are formed at the interface.

The usual damping effect due to the energy dissipation by viscosity and thermal conductivity exists irrespective of whether the mode is the convection or the gravity wave, but, for larger horizontal wave lengths, the effect of the boundary layer can become more important. The jump in the thermal conductivity in the boundary layer can give rise to overstability of the gravity wave in agreement with Souffrin and Spiegel (1967). The jump in the viscosity provides a self-catalytic action for the unstable flow if the viscosity is assumed to be the nonlinear turbulent viscosity due to the motion itself. The effect, however, is not strong enough to overcome the usual viscous damping.     

Boundary layer models of ocean floor spreading

Summary The equations of conservations of momentum and energy scaled with the characteristic values of the mantle indicate the presence of the upper boundary layer to produce the estimated rate of the ocean floor spreading by convection and the importance of the frictional heating. The depth of the upper boundary layer can be estimated from the balance of the viscous force with the horizontal pressure gradient at the sea floor. It is of the orders of 100 km and becomes deeper for the Pacific than for the Atlantic Ocean and also with frictional heating than without it. The frictional heating increases the surface heat flow of the heat conduction by ten to twenty percent for the Pacific Ocean but only by a few percent for the Atlantic Ocean. The similarity solutions are determined for the temperature and horizontal velocity in the upper boundary layer. These solutions are expressed in power series of the variabley x ^–n , wherex, y, andn are horizontal and vertical coordinates and numerical constant, respectively. Both temperature and horizontal velocity within the boundary layer are higher for the Pacific than for the Atlantic Ocean. When a larger viscosity is applied, it causes the increase of horizontal velocity below the surface because of the surface boundary conditions of the finite velocity and of vanishment of the velocity shear. The higher horizontal velocity generates higher temperature because it advects hotter material from the mid-ocean ridge site. The direct effect of frictional heating on the temperature distribution of the similarity solution is almost negligible, since the shear zone is deep and near the lower boundary of the upper boundary layer. In the similarity solution, the surface heat flow which is increased by the frictional heating is given as the boundary value. The effect of the frictional heating is important below the mid-ocean ridge.     

Studying ice sheet stability using the method of separation of variables

Abstract

Stability of two-dimensional stationary continental and marine ice sheets is studied using perturbations of ice sheet surface elevation and the margin position about a steady-state solution. Ice accumulation rate at the upper surface is specified as a function of elevation and span. Ice is considered as a Newtonian fluid. Linearisation and separation of variables yield a non-homogeneous eigenvalue problem. In case when the margin moves, a necessary condition for the existence of a solution is very restrictive on the functions of the bed profile and accumulation rate. The method of separation of variables is applicable when the margin is assumed to be stationary. Singularity of the perturbation at the margin is treated using the method of matched asymptotic expansions. Numerical experiments show that dependence of accumulation rate on elevation contributes strongly towards ice sheet instability. Bed slope, ice thickness at the grounding line of a marine ice sheet and equilibrium line inclination are the main parameters determining the ice sheet reaction to a surface perturbation.     

Intermittent convection

Abstract

A theoretical analysis of pseudo two-dimensional, finite-amplitude, thermal convection is made for an infinite Prandtl number fluid which is subjected to a constant heat flux out of the top boundary and insulated at the bottom. For large Rayleigh numbers the convective flow becomes intermittent and the system is characterized by the following cyclic process: the formation of a thermal boundary layer by diffusion, the instability of this layer when it becomes sufficiently thick, the destruction of the layer by the convective flow, the dying down of the convection, and the reforming of the thermal boundary layer by diffusion. The periodicity and the horizontal wave number of the intermittent convective flow are found to be independent of the depth of the fluid layer but depend on the rate of cooling and the properties of the fluid.     

The impact of slope length on the discharge of sediment by rain impact induced saltation and suspension

Simulations using a mechanistic model of raindrop driven erosion in rain‐impacted flow were performed with particles travelling by suspension, raindrop induced saltation and flow driven saltation. Results generated by both a high intensity storm, and a less intense one, indicate that, because of the effect of flow depth on the delivery of raindrop energy to the bed, there is a decline in sediment concentration, and hence soil loss per unit area, with slope length when particles are transported by raindrop induced saltation. However, that decline is reversed when the critical velocities that lead to flow driven saltation are episodically exceeded during an event. The simulations were performed on smooth surfaces and a single drop size but the general relationships are likely to apply for rain made up of a wide range of drop size. Although runoff is not always produced uniformly, as a general rule, flow velocities increase with slope length so that, typically, the distance particles travel before being discharged during an event increase with slope length. The effect of slope length on soil loss per unit area is often considered to vary with slope length to a power greater than zero and less that 1·0. The simulations show that effect of slope length on sediment discharge is highly dependent on the variations in runoff response resulting from variations in rainfall duration‐intensity‐infiltration conditions rather than plot length per se. Consequently, predicting soil loss per unit area using slope length with positive powers close to zero when sheet erosion occurs may not be as effective as commonly expected. Erosion by rain‐impacted flow is a complex process and that complexity needs to be considered when analysing the results of experiments associated with rain‐impacted flow under both natural and artificial conditions. Copyright © 2009 John Wiley & Sons, Ltd.     

Wave-induced setup of the mean surface over a sloping beach

A new theoretical approach for the wave-induced setup over a sloping beach is presented that takes into consideration the explicit variations of the surface waves due to bottom slope and viscosity. In this way, the wave forcing of the mean Lagrangian volume fluxes is calculated without assuming that the local depth is constant. The analysis is valid in the region outside the surf zone and is based on the shallow-water assumption. A novel approach for separating the viscous damping of the waves from the frictional damping of the mean flow is introduced, where the mean Eulerian velocity is applied in the bottom stress for the mean fluxes. In the case where the onshore Lagrangian mean transport is zero, a new formula is derived for the Eulerian mean free surface slope, in which the effects of bottom slope, viscous wave damping and frictional bottom drag on the mean flow are clearly identified. The analysis suggests that viscous damping of the waves and frictional dissipation of the Eulerian near-bed return flow could lead to setup outside the surf zone.     

A comparison of two physics-based numerical models for simulating surface water–groundwater interactions

Problems in hydrology and water management that involve both surface water and groundwater are best addressed with simulation models that can represent the interactions between these two flow regimes. In the current generation of coupled models, a variety of approaches is used to resolve surface–subsurface interactions and other key processes such as surface flow propagation. In this study we compare two physics-based numerical models that use a 3D Richards equation representation of subsurface flow. In one model, surface flow is represented by a fully 2D kinematic approximation to the Saint–Venant equations with a sheet flow conceptualization. In the second model, surface routing is performed via a quasi-2D diffusive formulation and surface runoff follows a rill flow conceptualization. The coupling between the land surface and the subsurface is handled via an explicit exchange term resolved by continuity principles in the first model (a fully-coupled approach) and by special treatment of atmospheric boundary conditions in the second (a sequential approach). Despite the significant differences in formulation between the two models, we found them to be in good agreement for the simulation experiments conducted. In these numerical tests, on a sloping plane and a tilted V-catchment, we examined saturation excess and infiltration excess runoff production under homogeneous and heterogeneous conditions, the dynamics of the return flow process, the differences in hydrologic response under rill flow and sheet flow parameterizations, and the effects of factors such as grid discretization, time step size, and slope angle. Low sensitivity to vertical discretization and time step size was found for the two models under saturation excess and homogeneous conditions. Larger sensitivity and differences in response were observed under infiltration excess and heterogeneous conditions, due to the different coupling approaches and spatial discretization schemes used in the two models. For these cases, the sensitivity to vertical and temporal resolution was greatest for processes such as reinfiltration and ponding, although the differences between the hydrographs of the two models decreased as mesh and step size were progressively refined. In return flow behavior, the models are in general agreement, with the largest discrepancies, during the recession phase, attributable to the different parameterizations of diffusion in the surface water propagation schemes. Our results also show that under equivalent parameterizations, the rill and sheet flow conceptualizations used in the two models produce very similar responses in terms of hydrograph shape and flow depth distribution.     

Ground surface temperature history in southern Canada: Temperatures at the base of the Laurentide ice sheet and during the Holocene

We use temperature profiles from 7 deep (≈ 2000 m) boreholes located in southern Canada to infer ground surface temperature histories (GSTH) during the Last Glacial Maximum (LGM) and the Holocene. Visual inspection of the heat flow and of the reduced temperature depth profiles reveals significant regional differences with some sites showing conspicuous signs of post glacial warming, and other indicating only very small changes in ground surface temperature. These differences are confirmed by the inversions of the temperature profiles. The most prominent variations in GST are found at the Sudbury, Ontario, sites where the present ground surface temperature is high. With the exception of Sept-Iles, Quebec, the other sites only show moderate or no variation in GST. For all the sites, except possibly Sept-Iles, temperatures at the base of the ice sheet during the LGM were at or slightly below the melting point of ice. Temperatures might have been lower, a few degrees below 0 °C, at Sept-Iles. These results are consistent with field observations and model predictions suggesting high velocity basal flows in the ice sheet above the studied regions. These new data on basal temperatures will provide better quantitative constraints on glacier flow dynamics. The inversions give a chronology for the retreat of the ice sheet comparable to other proxies. Inversion and direct modeling show that, following the ice retreat, there was a warm period between 2 and 5 ka with temperatures 1–2 K higher than present. The inversion yields a time for this episode 1–2 kyr more recent than that inferred by other proxies for the Holocene climate optimum (HCO).     

Internal waves and temperature fronts on slopes

Time series measurements from an array of temperature miniloggers in a line at constant depth along the sloping boundary of a lake are used to describe the ‘internal surf zone” where internal waves interact with the sloping boundary. More small positive temperature time derivatives are recorded than negative, but there are more large negative values than positive, giving the overall distribution of temperature time derivatives a small negative skewness. This is consistent with the internal wave dynamics; fronts form during the up-slope phase of the motion, bringing cold water up the slope, and the return flow may become unstable, leading to small advecting billows and weak warm fronts. The data are analysed to detect ‘events’, periods in which the temperature derivatives exceed a set threshold. The speed and distance travelled by ‘events’ are described. The motion along the slope may be a consequence of (a) instabilities advected by the flow (b) internal waves propagating along-slope or (c) internal waves approaching the slope from oblique directions. The propagation of several of the observed ‘events’ can only be explained by (c), evidence that the internal surf zone has some, but possibly not all, the characteristics of the conventional ‘surface wave’ surf zone, with waves steepening as they approach the slope at oblique angles.     

A Model Study of the Strong and Weak Wind,Stably Stratified Nocturnal Boundary Layer: Influence of Gentle Slopes

With the exception of intermittency and waves, a brief review of the observed and modeled mean structure of the nocturnal boundary layer (NBL) is presented. The effect of gentle slopes on strong and weak wind NBL was investigated here using a one-dimensional model, with a simple correction term to account for the slope effects, identical to the one used by Brost and Wyngaard (1978). The study indicates that the wind profiles, temperature profiles and surface layer turbulence characteristics are extremely sensitive to the imposed geostrophic wind when small slopes are present especially for light winds. This is due to the complex interaction between the buoyancy driven slope flow and the imposed geostrophic wind that in turn influence the shear generation of turbulence. Finally, the current issues in the modeling of weak wind boundary layer are discussed.     

The errors in surface runoff prediction by neglecting the relationship between infiltration rate and overland flow depth

Many simplifications are used in modeling surface runoff over a uniform slope. A very common simplification is to determine the infiltration rate independent of the overland flow depth and to combine it afterward with the kinematic-wave equation to determine the overland flow depth. Another simplication is to replace the spatially variable infiltration rates along the slope i(x, t) due to the water depth variations h(x,t) with an infiltration rate that is determined at a certain location along the slope. The aim of this study is to evaluate the errors induced by these simplications on predicted infiltration rates, overland flow depths, and total runoff volume. The error analysis is accomplished by comparing a simplified model with a model where the interaction between the overland flow depth and infiltration rate is counted. In this model, the infiltration rate is assumed to vary along the slope with the overland flow depth, even for homogeneous soil profiles. The kinematic-wave equation with interactive infiltration rate, calculated along the slopy by Richard's equation, are then solved by a finite difference scheme for a 100-m-long uniform slope. In the first error analysis, we study the effect of combining an ‘exact’ and ‘approximate’ one-dimensional infiltration rate with the kinematic-wave equation for three different soil surface roughness coefficients. The terms ‘exact’ and ‘approximate’ stand for the solution of Richard's equation with and without using the overland flow depth in the boundary condition, respectively. The simulations showed that higher infiltration rates and lower overland flow depths are obtained during the rising stage of the hydrograph when overland flow depth is used in the upper boundary condition of the one-dimensional Richard's equation. During the recession period, the simplified model predicts lower infiltration rates and higher overland flow depths. The absolute relative errors between the ‘exact’ and ‘approximate’ solutions are positively correlated to the overland flow depths which increase with the soil surface roughness coefficient. For this error analysis, the relative errors in surface runoff volume per unit slope width throughout the storm are much smaller than the relative errors in momentary overland flow depths and discharges due to the alternate signs of the deviations along the rising and falling stages. In the second error analysis, when the spatially variable infiltration rate along the slope i(x, t) is replaced in the kinematic-wave equation by i(t), calculated at the slope outlet, the overland flow depth is underestimated during the rising stage of the hydrograph and overestimated during the falling stage. The deviations during the rising stage are much smaller than the deviations during the falling stage, but they are of a longer duration. This occurs because the solution with i(x, t) recognizes that part of the slope becomes dry after rainfall stops, while overland flow still exists with i(t) determined at the slope outlet. As obtained for the first error analysis, the relative errors in surface runoff volume per unit slope width are also much smaller than the relative errors in momentary overland flow depths and discharges. The relation between the errors in overland flow depth and discharge to different mathematical simplifications enables to evaluate whether certain simplifications are justified or more computational efforts should be used.     

Simulated effect of a forest road on near‐surface hydrologic response and slope stability

Forest management practices often result in significant changes to hydrologic and geomorphic responses at or near the earth's surface. A well‐known, but not fully tested, hypothesis in hillslope hydrology[sol ]geomorphology is that a near‐surface permeability contrast, caused by the surface compaction associated with forest roads, can result in diverted subsurface flow paths that produce increased up‐slope pore pressures and slope failure. The forest road focused on in this study is located in a steep forested, zero‐order catchment within the H. J. Andrews Experimental Forest (Oregon). A three‐phase modelling effort was employed to test the aforementioned hypothesis: (i) two‐dimensional (vertical slice), steady‐state, heterogeneous, saturated subsurface flow simulations at the watershed scale for establishing the boundary conditions for the catchment‐scale boundary‐value problem in (ii); (ii) two‐dimensional (vertical slice), transient, heterogeneous, variably saturated subsurface flow simulations at the catchment scale for estimating near‐surface hydrologic response and pore pressure distributions; and (iii) slope stability analyses, using the infinite slope approach, driven by the pore pressure distributions simulated in (ii), for assessing the impact of the forest road. Both observed and hypothetical rainfall events are used to drive the catchment‐scale simulations. The results reported here support the hypothesis that a forest road can have an effect on slope stability. The permeability contrast associated with the forest road in this study led to a simulated altering of slope‐parallel subsurface flow with increased pore pressures up‐slope of the road and, for a large rainfall event, a slope failure prediction. Copyright © 2005 John Wiley & Sons, Ltd.     

Radially convergent groundwater flow in sloping terrain

Abstract

The groundwater flow equation governing the elevation (h) of the steady-state phreatic surface in a sloping aquifer fed by constant recharge over a bi-circular sector is rhh′ ? r ² Bh′ + Pr ² ? PR ² = 0, where r is the radial coordinate, P is a constant involving recharge and aquifer properties, and B is the slope of the aquifer—bedrock boundary. The derived flow equation describes radially convergent flow through a sloping aquifer that discharges to a water body of fixed head. One important simplification is that in which the width of the bi-circular sector is constant, and the draining land becomes a rectangular aquifer. The bi-circular sector and rectangular-strip groundwater flow problems are solved in terms of implicit equations. The solutions for the steady-state phreatic surfaces depend on the ratio of recharge to hydraulic conductivity, the slope of the aquifer-bedrock, and the downstream constant-head boundary. Computational examples illustrate the application of the solutions.     

An experimental investigation of the effects of thawed soil depth on rill flow velocity

Water flow velocity is an important hydraulic variable in hydrological and soil erosion models, and is greatly affected by freezing and thawing of the surface soil layer in cold high-altitude regions. The accurate measurement of rill flow velocity when impacted by the thawing process is critical to simulate runoff and sediment transport processes. In this study, an electrolyte tracer modelling method was used to measure rill flow velocity along a meadow soil slope at different thaw depths under simulated rainfall. Rill flow velocity was measured using four thawed soil depths (0, 1, 2 and 10 cm), four slope gradients (5°, 10°, 15° and 20°) and four rainfall intensities (30, 60, 90 and 120 mm·h⁻¹). The results showed that the increase in thawed soil depth caused a decrease in rill flow velocity, whereby the rate of this decrease was also diminishing. Whilst the rill flow velocity was positively correlated with slope gradient and rainfall intensity, the response of rill flow velocity to these influencing factors varied with thawed soil depth. The mechanism by which thawed soil depth influenced rill flow velocity was attributed to the consumption of runoff energy, slope surface roughness, and the headcut effect. Rill flow velocity was modelled by thawed soil depth, slope gradient and rainfall intensity using an empirical function. This function predicted values that were in good agreement with the measured data. These results provide the foundation for a better understanding of the effect of thawed soil depth on slope hydrology, erosion and the parameterization scheme for hydrological and soil erosion models.     

Flow of a double-diffusive stratified fluid in a differentially-rotating cylinder

Abstract

A study has been made of a basic state of axisymmetric flow, at large rotational Reynolds numbers, in a double-diffusive stratified fluid contained in a vertically-mounted, differentially-rotating cylindrical cavity. The aim is to describe the qualitative characteristics of the flow of a fluid, the density of which is stratified by two diffusive effects, i.e., temperature and salinity gradients. Attention is confined to situations in which the temperature and salinity gradients make opposing contributions to the overall density profile, the undisturbed stratification being gravitationally stable. Finite difference numerical solutions of the governing Navier-Stokes equations have been obtained using the Boussinesq approximation. The results are presented in a way that illustrates the explicit effects of double-diffusivity when the cavity aspect ratio, height/radius, is O(1). The principal non-dimensional parameters characterizing the flow field are identified. In the interior core, the primary dynamic balance is between the horizontal density gradient and the vertical shear of the prevailing azimuthal velocity. The effective stratification is seen to decrease as the double-diffusivity increases, even if the overall stratification parameter, St, is held constant. The solute field contains a very thin boundary layer structure at large Lewis numbers. The effective stratification increases with the Prandtl number. Results have been derived for extreme values of the cavity aspect ratio. For small cavity aspect ratios, the dominant dynamic ingredients are viscous diffusion and rotation. For large aspect ratios, the bulk of the flow field is determined by the rotating sidewall. In this case, the direct influence of the double-diffusivity is minor.     

A nonlinear,resistive boundary layer in rotating hydromagnetic flow

A new nonlinear boundary layer in rotating hydromagnetic flows is presented. The purpose of this layer is to provide a smooth transition between the magnetic field at a rigid, electrically insulating boundary and that far from the boundary. The boundary layer problem is solved in a cylindrical geometry, assuming that the variables obey the hydromagnetic analog of the von Karman similarity. The primary momentum balance within this boundary layer is between pressure, Coriolis and hydromagnetic forces; the viscous and nonlinear inertia forces are unimportant. Diffusion of the magnetic field is balanced by advection where the advecting fluid flow is induced by the action of the hydromagnetic body forces within the boundary layer. The structure of the layer depends upon a single parameter which measures the strength of the normal magnetic field. Steady solutions are presented and their uniqueness and temporal stability are analyzed. The relevance of this layer to the core of the Earth is discussed. It is estimated that this boundary layer would be at least as thick as one tenth of the core radius if it exists within the core.     

Array measurements of the bottom boundary layer and the internal wave field on the continental slope

Abstract

An array of current meters was placed on the continental slope and rise for two months in the autumn of 1970. The bottom boundary layer was penetrated on the slope. On the smallest array scale, of the order of 1 kilometer, the array functioned as a directional internal wave antenna. Moving shoreward, the current spectra show strong suppression of the inertial peak and strong enhancement of the semidiurnal tide. The measured wave number spectra show that the tidal energy is almost completely baroclinic, and probably being generated in the region where the slope becomes “critical” for the tidal period. If this area is typical of worldwide conditions, a substantial fraction of the dissipation of surface tides takes place on the continental slopes by conversion to baroclinic waves. The bottom boundary layer has been modeled by an extension of the work of Ellison (1956) to a sloping boundary in a fluid of positive stability. An equivalent constant eddy coefficient has the value 3 cm²/sec as determined from the measurements.