Introduction
In 1957, the Journal of Glaciology published the iconic paper that set an important scientific direction for glaciological research that endures to this day. Johannes Weertman's ‘On the sliding of glaciers’ was, according to Fowler (Reference Fowler2010),
[t]he first paper in the subject which properly addresses the physics….
This paper (Weertman, Reference Weertman1957a) ignited follow-up research by many of the world's most prolific glaciologists over the subsequent six decades, and is one of the most widely cited in the history of glaciology (over 640 citations listed in Google Scholar, more than 50 in 2018–19).
The theory of glacier sliding, as envisioned by Weertman in 1957, endures as the pedagogical source for understanding the phenomenology of ice sheet and glacier sliding even though, it involves an unconventional representation of bedrock roughness, and only two of the five physical processes that are involved. The theory accounted for enhanced creep motion and flow by regelation, but did not account for bed separation/cavity formation, hydraulic pressure effects and deformation of the material below the ice bed interface. Weertman's 1957 glacier sliding problem has been revisited several times before (e.g. as described in Fowler, Reference Fowler2010). The study presented here will not be as comprehensive, or as critical, as some (e.g. Lliboutry, Reference Lliboutry1968), because my intention is to merely provoke further research in the problem. My aim is to revisit Weertman's 1957 glacier sliding problem as originally posed with its distinctive tombstone bedrock obstacles using numerical methods. The goal is to both celebrate the achievement embodied in his paper and to identify several elements of the problem that may warrant further research.
Weertman's early life and origin of his interest in glaciology
Johannes ‘Hans’ Weertman (Fig. 1) grew up in Beaver, Pennsylvania (USA), as the oldest of two sons of Roelof and Christina Weertman, who emigrated to the USA from Rotterdam, The Netherlands, in 1923, 2 years before Hans's birth. When Hans was growing up, Roelof was an engineer and mechanical draughtsman working for the St. Joseph Lead Company of Pittsburgh, PA. In his spare time, and later in retirement, Roelof made violins and cellos. During Hans's high school days, Roelof built a violin (Fig. 2) that was taken by a Beaver, PA neighbor just a few years older than Hans, Harrison Richardson, to the Ross Ice Shelf, Antarctica, as part of the third and final American Expedition led by Richard Byrd. The violin wintered over at Little America III and undoubtedly provided the expedition with artistic expression. This was the earliest intersection between Hans and the glaciological world where he was to make key scientific contributions in his later life, i.e. notably including his other 1957 paper on the flow of ice shelves (Weertman, Reference Weertman1957b) and later papers on stability of the West Antarctic Ice Sheet (the first being Weertman, Reference Weertman1974).
Graduating from high school in 1942, Hans joined the US Marines as a private to serve in World War Two. Hans was assigned to a Navy V-12 unit, an officer's training school, at Pennsylvania State University from 1942 to 1945 before serving in the Pacific, where he was deployed as a corporal in the 2nd Marine Division in summer of 1945. He spent from September 1945, to January 1946, with the 2nd Division Headquarters battalion in Nagasaki, Japan. Following his discharge in spring of 1946, he completed his education leading to a doctorate of science in 1951 from Carnegie Institute of Technology (now Carnegie Mellon University).
Hans' graduate study was supervised by James Koehler, a notable solid-state physicist and materials scientist who was among the first in the USA to embrace the dislocation theory for the permanent deformation of solids. There was no Earth science or glaciology in Hans' education, as Carnegie Institute of Technology had no Earth-science department. This is why the first 17 papers he published (listed in his CV provided to me by Northwestern University) were either on the phenomenology of specific metallurgical creep examples (9 papers) or on the burgeoning theory of dislocation creation and movement in metals (eight papers).
Casual examination of the figures and equations displayed in these initial 17 papers gives only one clue to Weertman's burgeoning interest in ice. The fourth paper he published, ‘Theory of steady-state creep based on dislocation climb’ (Weertman, Reference Weertman1955), was received and accepted in 1955, and was the first to reference anything glaciological. It cited John Glen's experimental results (Glen, Reference Glen1955) on creep deformation of ice. Weertman's citation of Glen's work was in the last entry of a table displaying creep-test data on metals, including aluminum, lead, tin, platinum, bromine and silver.
Weertman developed the tombstone bed sliding model in 1955 and 1956 while he and his wife, Julia Randall Weertman, were employed by the U.S. Naval Research Laboratory (NRL) in Washington, DC. In 1982, he spoke of his experience with developing the tombstone bed model (Weertman, Reference Weertman1982):
I worked on it mainly on weekends while driving on the Pennsylvania Turnpike to our parents' homes near Pittsburgh. My wife knew when I was thinking about it because of an ever slowing down car, a sign for her to take over driving.
He went on to say:
The story goes on. Peter Haasen visited NRL and told me that John Glen earlier had found power law creep in ice. This led me into the glaciological literature and I learned that how glaciers slide was not understood. From my thermo course, I dimly remembered reading in Fermi's book that glaciers slide around bed bumps by the pressure melting phenomenon. A simple calculation showed me that this is too slow a process for large bumps. But it was immediately clear that creep lets ice get around big bed bumps quickly. So I was able to write the first theory on glacier sliding.
Peter Haasen was a scientist at the Metallurgical Laboratory of the University of Chicago at the time, and went on to become a distinguished professor of metal physics at the Georg August University of Göttingen. Since his interaction with Weertman, however, he is cited only obliquely in the glaciological literature making reference to his work on recrystallization and fabric development in metals.
Weertman's glacier-sliding problem
When I first saw Figure 1 of Weertman's first paper on a glaciological subject (replicated as Fig. 3 here), I was struck by his unconventional representation of the glacial bed geometry, which convention tends to portray as smoothly ground down by ice flow. Only now do I suspect that the tombstone obstacles Weertman envisaged were salted with genius, because that specific geometry swept aside, in a single stroke, many, if not all, of the mathematical inconveniences that bar a simple, intuitive understanding of how glaciers flow across their bed. Weertman was well aware of the fact that glacier beds were eroded and smooth at the time of his paper, because he had visited the valleys and glaciers of Switzerland in 1951/52 while a Fulbright post-doctoral fellow at Ecole Normal Supérieure, Paris. He had even walked into the U-shaped Furka Pass in 1955, as his family photos attest. He chose to adopt the sharp-cornered tombstone because he believed that this idealization would allow him to reach qualitative and quantitative conclusions that were universal. He may have further been drawn toward the rectilinear geometry because of his background in metallurgy and solid-state physics, where fixed angles of crystal structure defined both rock and ice on a microscopic scale.
Once Weertman postulated the unconventional bed geometry, he made quick work of arriving at a useful working solution of the glacier-sliding problem. He defined it as
$\dots$ simply the determination of the velocity U of this block of ice (Weertman, Reference Weertman1979)
while referring to Figure 1 of his 1979 paper which depicted the force balance of a block of ice being dragged by the glacier above across a bedrock surface. In the 1957 paper, Weertman assumed that an average horizontal drag force τl 2 was experienced by the ice in each unit cell of the ice/bed representation. Here, l is the horizontal dimension of the unit cell (tombstone spacing) and τ is a basal shear stress experienced at large scales by the glacier above. Gravity was not accounted for in the 1957 paper, as its importance for the ice near the bed was dismissed in relation to the stress applied from the ice above. He assumed that the bedrock surfaces were perfectly lubricated by meltwater, and thus could not sustain tangential forces. This assumption, applied to the tombstone geometry, determined the difference, τl 2/a 2, between normal forces on the stoss and lee faces of each tombstone. Here, a is the dimension of the tombstone obstacle. This force difference gave the two key ingredients of Weertman's theory: a means to determine (1) temperature difference, heat flux and thus the meltwater flow around the obstacle (by regelation flow), and (2) the degree to which stress, and therefore strain rate, would be enhanced in the vicinity of the obstacle.
Weertman's 1957 understanding of regelation sliding was inspired by a brief passage in Enrico Fermi's textbook on thermodynamics in a section describing the Clapeyron equation:
When the mass of ice encounters a rock on the glacier bed, the high pressure of the ice against the rock lowers the melting point of the ice at that point, causing the ice to melt on one side of the rock. It refreezes again immediately after the pressure is removed. In this way the mass of ice is able to flow very slowly around obstacles (Fermi, Reference Fermi1937, chapter 4, section 15).
Weertman was immediately wary that regelation sliding alone would not explain the high rates of sliding observed, because he quickly surmised that its effect would be negligible when bedrock obstacles exceeded the centimeter scale. He nevertheless proposed that the pressure difference between the leading and trailing faces of the tombstone obstacle would, by virtue of the co-existence of ice and water at the bed, and the effect of pressure on melting temperature, produce a heat flux from the downstream side to the upstream side of the obstacle within the bedrock. He did not consider the full thermodynamic problem of how the bedrock surface would be maintained at the pressure melting temperature, and simply assumed this condition to be true. He further simplified the reasoning necessary to deduce what specific pressure would apply in the determination of the melting temperature by conducting a simple force balance, assuming pressure to be constant on leading and trailing faces of the obstacles. The thermodynamically complete statement of the regelation problem came in the decade after Weertman's 1957 study, most notably in the work of Nye (Reference Nye1967, Reference Nye1969, Reference Nye1970) and Kamb (Reference Kamb1970).
Weertman then reasoned that, since a < l, the stress field in the ice in the immediate vicinity of each tombstone should be amplified over the ambient scale τ that applied to the ice above the basal layer. This amplification, combined with the > 1 exponent in the flow law for ice, meant that basal ice would creep around each tombstone at a rate that was enhanced relative to the general background creep of the ice above the bed. Weertman calculated the amplified strain rate near each tombstone, and then reasoned:
To obtain a velocity of sliding from the creep rate it is necessary to know the distance over which the stress is essentially at the level ${\displaystyle {1} \over {2}}\tau l^2/a^2$. A reasonable estimate which uses Nye's (Nye, Reference Nye1953) values of the stress around contracting holes gives the value a, on either side of the obstacle, for this characteristic distance.
This then gave him the estimate of the glacier sliding flow due to creep. Weertman (Reference Weertman1957a) later pointed out:
The sliding velocity caused by the stress concentration [surrounding the obstacles] is equal to the [enhanced] creep rate times the distance in the direction of motion over which this creep rate is effective. The sliding velocity is, of course, only a rough estimate [because] we do not know exactly the distance over which the [enhanced] creep rate is effective.
This probably reflected his recognition that his simplified consideration of a basal layer of arbitrary dimensions would not ultimately provide a quantitatively accurate depiction of glacier velocity. At the time he wrote this, relatively little attention was paid to how a sliding velocity at the bed would combine with deformational flow above the bed to yield the surface velocity of a glacier. This was possibly because there were few observations of glacier flow before 1961 relative to today (e.g. from remote sensing methods such as by Rignot and others, Reference Rignot, Mouginot and Scheuchl B2011; Mouginot and others, Reference Mouginot, Scheuchl and Rignot2012).
A numerical solution of the problem
In the current study, the velocity of ice in the vicinity of the glacier bed is partitioned into creep and regelation components in the same manner as the original study (Weertman, Reference Weertman1957a):
The creep component of the sliding velocity, uc, is determined using the finite-element method to solve the steady-state stress-balance equation for an incompressible, Newtonian viscous rheology. I defer using the power-law rheology used by Weertman (Reference Weertman1957a) to future study, because results obtained with a Newtonian rheology appear to be sufficiently interesting for the purpose of the current study. Assuming small Reynolds number and excluding gravity (as Weertman assumed was justified for characterizing the velocity in the immediate vicinity of the bed), the governing equation for uc is:
where p is the pressure perturbation due tothe effects of sliding near the bed, assumed to be additive to glaciostatic pressure from the ice above, and η is the viscosity, assumed constant. For boundary conditions, free-slip with zero-tangential stress is assumed on all bedrock surfaces and a spatially uniform shear stress τxz = τ is applied to the top of the numerical domain, which is taken to be the z = l surface, where z = 0 constitutes the plane on which the tombstone obstacles are arrayed (Fig. 3), and where x and y are the horizontal coordinates. It is noted that the disregard of gravity (i.e. assuming g = 0) requires that the thickness of the numerical domain, here taken to be l, be sufficiently small to avoid complexity associated with changing stress over the vertical extent of the numerical domain. In the current study, the vertical extent of the numerical domain could have been taken to be somewhat larger than the horizontal spacing of tombstones (l) and still have satisfied the criterion of being sufficiently small. The choice made to limit the vertical extent to a distance l above the bed is also made in this study to limit the number of nodes in the finite-element mesh, thereby reducing the computational burden. Had z = 2l been chosen, for example, the results would be quantitatively different, as this alternative choice would imply that τxz = τ is applied as a boundary condition at a greater distance from the obstacles. Simulations conducted (not shown) to examine the sensitivity of results to the vertical extent of the numerical domain suggested that qualitative results, which are the focus of this study, are adequately portrayed when l is the vertical dimension of the domain.
As implied by the regular geometry, symmetry conditions apply on vertical planes that bound the numerical domain horizontally, and which cut the unit cell of the tombstone array as shown by the red line in Figure 3. Dimensionless variables are used:
with the relation
to simplify the governing equation (the viscosity parameter disappears) and to render the results of the numerical solution independent of the chosen Newtonian viscosity parameter η. For the simulations shown in this study the aspect ratio of the tombstone obstacles, ν = a/l was taken to be 1/4.
Results
The solution, uc was evaluated using the finite-element simulation software package COMSOL, and is displayed in Figure 4. What was not anticipated by Weertman in his 1957 paper was the complexity of flow immediately in front of the upstream (and by symmetry, downstream) face of the tombstone obstacle. The flow complexity is also illustrated in Figure 5, where several specific streamlines that impinge on the upstream face of the tombstone obstacle undergo inversion (with respect to z). An equivalent, but opposite, deformation of the streamlines occurs on the downstream side of the obstacle.
The illustrated flowline geometry (Fig. 5) for tombstone obstacles suggest that striations on the leading edge and trailing edges of the obstacle, and in its immediate vicinity on the bed will deviate considerably from the general direction of flow in the ice above the bed. This is something that should be kept in mind when interpreting striation direction on formerly glaciated terrains. For an accurate interpretation of striation direction as an indicator of large-scale flow direction, striation direction should be measured in an area that is large compared to bed obstacle size, and this thus necessitates awareness of the geometry of bedrock obstacles in the field area. Simulations conducted using smooth, hemispherical-shaped bedrock obstacles (one of which is shown in the Appendix) revealed that the complexity of the implied striations (based on flowline directions) is reduced when bedrock lacks the sharp, angular edges of tombstones. This means that formerly glaciated terrains that are smooth due to a long period of glacial erosion are less likely to display striations that deviate from the large-scale direction of former ice flow.
Had Weertman been more interested in the fluid-dynamic details of the flow around the tombstone obstacle he may have encountered the work of Moffatt (Reference Moffatt1964) which suggests that recirculating and twisting flows may be quite common on inside corners of obstacles imbedded in a viscous flow. Meyer and Creyts (Reference Meyer and Creyts2017) vividly illustrate a compelling glaciological application of Moffatt's work in the case of flow within sub-glacial valleys in the Gamburtsev Mountains of Antarctica (see also Gudmundsson, Reference Gudmundsson1997). The mere presence of the twisting flow upstream and downstream of bed obstacles may be of minor importance in the development that Weertman was trying to accomplish; however, such considerations may be of great interest in ice-core problems where preservation of stratigraphic order is important.
Another aspect of creep-deformation sliding not discussed in Weertman's 1957 paper, is that both pressure p and vorticity $\vec \omega = \nabla \times {\bf u}_{\rm c}$ appear to have an unphysical singularity along leading and trailing edges of the tombstone obstacle. On the upstream side of the obstacle (Fig. 6a), a curious ‘double stripe’ pattern is displayed for pressure p. Red color (high positive pressure) lines the edge on the face which opposes the flow, and blue color (high negative pressure) lines the edge of the faces which form the top and side of the obstacle. While not shown, the downstream edges that frame the downstream face of the tombstone show a pattern of opposite sign. No matter how the finite-element mesh was refined to increase spatial resolution in the neighborhood of the edges, the stripe-like pattern of pressure appeared at the minimum resolvable length scale, shrinking ever more closely to the 90° corners of the tombstone edges. Singularities of a similar nature, but caused by effects associated with cavity formation over smoother, sinusoidal bedrock obstacles were seen, for example, in the numerical study by Gagliardini and others (Reference Gagliardini, Cohen, Raback and Zwinger2007). It is noted that total pressure in ice cannot be arbitrarily large and cannot be negative without phase changes (for extremely large positive pressure) or bed separation or tensile fracture (for even modest negative pressure).
The cause of the singularity is the infinite curvature of the sharp, 90° edges of the tombstone. On these edges, vorticity magnitude $\vert \nabla \times {\bf u}_{\rm c} \vert$ must be infinite if non-zero tangential flow exists on each bedrock face leading to the edge. That the pressure should have a singularity when the vorticity does is a well-known property of viscous flows (e.g. see Barcilon and MacAyeal, Reference Barcilon and MacAyeal1993). The presence of the singularity in p has important consequences on the computation of ur.
I computed ur by simulating the heat transport in the bedrock assuming that the bedrock surface is everywhere at the pressure melting temperature. I use the heat-conduction equation in the domain below the bed bounded by z = −l (chosen, again, for computational reasons) at the bottom, and by the ice/bedrock interface at the top (Fig. 3). The effects of geothermal heat flux, atmospheric temperature at a distant upper surface of the glacier, strain-heating and the presence of an englacial meltwater component are disregarded. In their place, I assume the numerical domain to be surrounded by ice at the pressure melting temperature.
The heat-conduction equation with constant thermal conductivity was solved with boundary conditions formulated to account for symmetry (y = 0, l/4), antisymmetry (x = 0 − l/4), and thermal insulation at z = −l (again, a computationally motivated simplification). At the ice/bedrock interface, temperature T was specified using the Clapeyron equation:
where S is the locus of the bedrock surface, T m denotes the pressure-dependent melting temperature, T o = 273.15 K, C = 7.4 × 10−8 kPa−1 is the pressure-dependent melting temperature, and p is the pressure perturbation derived from the solution of the creep sliding problem discussed above, i.e. associated with uc. Weertman (Reference Weertman1957a) did not solve the stress balance equations to deduce a pressure variation on the leading and trailing faces of the obstacle, but rather deduced them from a simplified consideration of total force balance. Here, the pressure is taken from a solution of the stress balance equations used to determine uc. In this way, my study differs from Weertman's.
The heat flux associated with the steady state solution of the heat-conduction equation is shown in Figure 7. As anticipated by Weertman, the computed heat flux normal to the ice/bedrock interface on the upstream side of the tombstone is indeed positive, and by symmetry, negative on the downstream side of the tombstone. The normal heat flux, integrated over the leading and trailing faces of the obstacle were approximately the same as deduced by Weertman (differing only in the fact that his estimate of pressure was simplified relative to the numerical solution), and the conclusion he reached concerning the inverse dependence of this heat flux on obstacle size (i.e. depending on 1/a) was confirmed by the numerical solution.
The most important difference between the numerical solution and Weertman's simplified heat flux, however, is due to the singularity in p associated with the creep sliding flow (uc). This singularity in the numerical solution gives a strongly non-uniform pattern of heat flux normal to the boundaries of the tombstone. In addition, there are strong normal heat fluxes on the top and sides of the tombstone which are not perpendicular to the direction of flow. Another difference is that there are broad regions of relatively weak normal heat flux on the flat horizontal parts of the bedrock interface between the tombstone obstacles. If, as an alternative interpretation, the normal traction to the bedrock surface n · (T · n), where n is the unit normal to the bedrock surface and T is total stress, is used instead of p in Eqn (7), the non-uniform pattern of heat flux on the tombstone boundaries remains qualitatively the same. Thus, a question of whether the pressure or total normal stress applies in the Clapeyron equation does not matter in determining this result.
If the computed heat fluxes are used to estimate the boundary condition applying to ur at the ice/bedrock interface, the resulting regelation flow is not spatially uniform, and is thus not a ‘plug flow’. This means the creep flow must be coupled to the regelation flow. I did not attempt the challenging numerical implementation of such a coupling in the current study.
Conclusion
Having worked solely on metallurgical problems until 1956, Weertman (Reference Weertman1957a) approached the glacier-sliding problem as a novel Earth-science curiosity that exemplified the burgeoning field of dislocation mechanics in solids. In the fourth sentence of the introduction to his 1957 paper, he invokes the commonality of metallurgy and glaciology:
The creep behavior of metals is quite similar to that of ice. For this reason the extensive work on the creep of metals is of interest to glaciologists and the work on glaciers is interesting to metallurgists.
The progress Weertman made toward understanding glacier sliding in terms of dislocation mechanics was remarkable, but Weertman's work did not go so far as to present a complete embodiment of the three-dimensional flow field. Scientists who did embrace this three-dimensional structure were forced to depart from Weertman's tombstone problem and consider two-dimensional flow fields with smooth, continuously differentiable bed geometries that were amenable to mathematical analysis (e.g. as recounted by Fowler, Reference Fowler2010) and more conforming to what is expected when glacial erosion is active. In the current study, I have provided a partial view of the three-dimensional flow field in the tombstone problem that are only addressable with computational methods. Indeed, computational methods at the time of Weertman's study were out of the reach of glaciology considering that the application of computational methods to geophysical problems had only begun in 1950 (Charney and others, Reference Charney, Fjörtoft and vonNeumann1950). Finally, it is notable that one of the first computer simulations used in glaciological research was by Almut Iken (Reference Iken1981), the purpose of which was to investigate bed separation and cavity formation in a bed geometry which also had sharp corners.
Three aspects of Weertman's glacial sliding problem are revealed in the numerical solutions presented here: (1) the existence of stratigraphy-inverting flows in the interior corner regions just upstream and downstream of the tombstone obstacles, (2) pressure and vorticity singularity on the edges of the tombstone, and (3) the failure of the presumption that regelation flow around sharp-cornered obstacles can be plug-like, and thus uncoupled from creep flow. Whether it is necessary to understand these three aspects remains to be determined, as the tombstone-bed model of glacier sliding remains an unconventional idealized stepping stone needed for pedagogical understanding. It is, however, concluded that for both students and advanced researchers, revisiting the development of Weertman's 1957 theory, or of any classical theories in glaciology, provokes further questions.
Acknowledgments
This work was supported by U.S. National Science Foundation grants PLR-1443126 and NSFGEO-NERC1841467. I thank Bruce Weertman for sharing information about his father's life and glaciological work and for providing the two family photographs that appear herein. Comments by Tim Creyts, Colin Meyer, and Alan Rempel were helpful with early stages of the manuscript. Two anonymous referees and Scientific Editor Martin Truffer corrected several misunderstandings I had, and suggested citations to additional aspects of the history of the problem, notably the advances by Nye and Kamb in the years following the 1957 study. I especially thank the Scientific Editor for discussing the concepts involved in the manuscript and for conducting thorough proofreading. I thank the Northwestern University Department of Materials Science and Engineering for inviting me to the symposium and remembrance celebrating the lives and contributions of Johannes and Julia Weertman where I gathered personal anecdotes attributed to Johannes, and where I met other members of his family.
APPENDIX A. Glacier sliding associated with hemispherical obstacles
To provide a contrast to the results for tombstone obstacles, simulation of uc and ur using hemispherical obstacles were conducted. The hemispherical obstacles have radius $r = \sqrt {{2}/{\pi }}a$ giving an equal cross-sectional area (a 2) as the tombstone. The large-scale pattern of uc for the hemispherical obstacle (Figs. 8a and b) was found to be similar to that of the tombstone obstacle (Fig. 4). However, there were no singularities in p (as is evident from the pressure-determined melting temperature shown in Fig. 9a) or vorticity (not shown). Streamlines (Fig. 8c) did not show the complex overturning characteristic of the tombstone obstacles (Fig. 5).
In the computation of ur for the hemispherical obstacle, as for the tombstone obstacle (Fig. 7b), non-zero heat flux normal to the bed (not shown) was found on parts of the upward-facing boundary between the hemispherical obstacles when the perturbation pressure p was used to determine the melting temperature using Eqn (7). This means plug-like regelation flow, which is everywhere parallel to the z = 0 plane, is not possible even for hemispheres. Melting and freezing on parts of the bed between the hemispherical obstacles generate the need to couple the regelation process with the creep process and refute the presumption of a plug-like flow.
To further investigate the enigmatic absence of a plug-like regelation flow even for bedrock obstacles that prohibited singularities in pressure, I conducted an additional simulation where the melting temperature was prescribed according to the suggestion by Nye (Reference Nye1967):
where T is the total stress tensor and n is the outward pointing normal to the ice/bedrock interface. This modification accounts for the fact that pressure in the water film is different from that in the ice, and must be in balance with the total normal stress at the contact between ice and water film (see Rempel and Meyer, Reference Rempel and Meyer2019, for further discussion).
Using Eqn (A1) substantially changed the simulated temperature field on the bedrock interface (Figs. 9b and c), and removed the bed-normal heat flux on flat parts of the bed between the obstacles. The simulated bed-normal heat fluxes were somewhat marred by numerical noise (which is less visible in the tombstone geometry simulation owing to the dominance of the singularity at the corners), but had the large-scale pattern that is more consistent with plug-like regelation flow.
As shown in Figure 10a, the heat flux normal to the ice/bedrock interface is mostly restricted to the spherical obstacle, and mimics the pattern of the value of the x-component of the normal vector to the ice/bedrock interface (Fig. 10b). Heat flux trajectories associated with the improved simulation are shown in Figure 10c. The heat flux pattern in the bedrock is more like what Weertman envisioned for the tombstone obstacle, simply connecting the upstream with the downstream faces of the obstacle. The recirculating heat flux trajectories around the edges where the hemisphere connects with the basal plane at z = 0 suggest that the numerical solution does not support an exactly plug-like regelation flow. Whether this is a universal result or simply an artifact associated with finite numerical resolution is relegated to future study.