Energy balance model

Our energy balance model is designed to capture the essential physics underpinning the sinking process of an englacial meteorite while avoiding extraneous details. The model does not attempt to offer highly accurate Antarctic site-specific results; it aims to address the proof-of-concept question of englacial sinking through the inclusion of dominant (or aggregated) terms in the various energy balance equations and so should be taken to offer generalized quantitative predictive results at this stage.

We consider three distinct stages within our modelling. Chronologically, the evolution starts with the meteorite heating due to absorption of downwelling solar radiation through the ice; however, it remains fully encased in ice while temperatures remain below freezing throughout. Then, once the top surface of the meteorite is sufficiently heated it is able, through conduction, to induce melting of the overlying ice, allowing an upper water layer to form; this water layer can itself evolve through its interaction with both the overlying ice and the underlying meteorite. Next, if the meteorite continues to warm sufficiently, melting at its lower surface will also commence. When it does, we allow the sub-meteorite melt water to be squeezed upwards by the weight of the meteorite adding to the overlying water layer, thereby enabling the meteorite to sink downwards and fill the displaced volume. These three stages can stop at any time, and freezing reoccur, as the meteorological inputs vary, thus bringing the sinking process to a halt (this being the case for the Antarctic situation during the winter months).

To model this process in the simplest and most elucidating manner, we consider a one-dimensional representation to the physical problem, as sketched in Fig. 3. Thus, the model has four distinct regions at increasing depth, z: (i) an upper ice layer that is exposed to the atmosphere; (ii) a water layer (which need not always exist); (iii) the meteorite; and (iv) a lower ice layer. To model the Antarctic situation, one must incorporate the relative motion of the meteorite to the ice surface. With the MSZ surface in equilibrium at z=0, say, the upwelling ice velocity V must be matched by the ablation rate (which is the sum of the energetic sublimation rate v, and the non-energetic rate that ice is scoured off the surface by the katabatic winds, V−v).

Figure 3: Model geometry. A (not to scale) schematic diagram highlighting the boundaries and geometry of the mathematical model for the Antarctic situation, in which an englacial meteorite is exposed to solar radiation. Full size image

To frame the problem mathematically, we need to consider the energy balance at the five boundaries/interfaces of the four regions. It is of note that two of these interfaces are ‘free-boundaries’, whose locations are variables that require solving as part of the problem: the upper ice/water interface (z=a), and the water/meteorite interface (z=b). The time-varying solution to these two variables thus determines the dynamics of the meteorite within the ice. In solving for them, we must also solve for the other variables involved within the model, namely the temperatures T j in regions j=1−4 (note that the air temperature is a model input).

At the atmosphere/ice interface (z=0), the energy balance is given by16:

where, sequentially, the terms on the left hand side represent: the contribution of the shortwave solar flux; the incoming longwave flux; the linearized outgoing longwave radiation (Stefan–Boltzmann’s law); the sensible heat flux; and sublimation, respectively. The right hand side is the heat flux into the ice, which is given by Fourier’s heat-transfer law. All model parameters are defined in Table 2.

At the ice/water interface (z=a; once it exists) the temperature must be at 0 °C. Further, the energy balance tells us that the energy flux for melting the ice (or its refreezing) plus the heat flux into the ice must equate with the heat flux from the water layer:

where the overdot denotes differentiation with respect to time. Note that the inclusion of +V on the left hand side is due to the frame being fixed with respect to the upwelling ice. When the water layer does not exist, we neglect the phase change term in (2) and k w becomes k i .

At the water/meteorite interface (z=b), we must consider the shortwave solar energy that is reaching the meteorite surface (the longwave flux having been absorbed at the surface), and how this flux decays with depth below the surface. To achieve this we make use of the Beer–Lambert’s exponential decay law with extinction coefficient γ i . The balance at this interface is given by the solar energy flux plus the heat flux from water layer equating with the heat flux into the meteorite:

Here Q r may be thought of as the shortwave energy that would be absorbed by the meteorite in the absence of absorption/scattering by the ice, and thus Q r =(1−α i )(1−α m )(1−s h )S, where S is the total incoming solar flux, α i is the ice surface broadband albedo, α m is the broadband meteorite albedo and s h is the percentage of incoming solar radiation that is shaded (thus capturing the effects of local topography). Due to the inclination of the sun to the horizontal (solar elevation angle θ), we make use of an effective distance b′ that represents the distance travelled while being attenuated within the ice and melt water (where we neglect the tiny effect of refraction within the very-thin melt-water layer):

At the lower meteorite/ice interface (z=b+w), where w is the meteorite width, one needs to compute whether the temperature is high enough to melt the ice. When the interface temperature is at 0 °C, melting is permitted, and the energy balance is given by the heat flux from the meteorite layer equating with the energy flux for melting plus the heat flux into the lower ice layer:

Again, the +V term is added due to the moving frame of reference. When the temperature is below melting point, we neglect the phase change term in (5).

At the bottom of the lower ice region we shall prescribe a temperature-gradient condition of the form

In the laboratory, where the ice block is relatively warm and thin (5 cm), and the incoming heat fluxes held constant, we use a measured value of the ice temperature T ∞ at z ∞ . Thus, the ‘far-field’ temperature gradient is

However, in the Antarctic setting we are unable to prescribe such a measured temperature at z ∞ , as the temperature profile within the ice varies with time (albeit slowly). As an alternative, we note that the ice underneath the meteorite must tend to the temperature profile of the surrounding ice, and so we must prescribe a condition drawn from the meteorite-free situation, that is, we need to match the temprature flux at z ∞ with that of the ice thermal boundary layer. To achieve this, we can utilize the fact that the magnitude of the (meteorite-free) annual ice-temperature variation in the boundary layer diminishes rapidly with depth, which means that in practice there is only a small average annual temperature gradient at z ∞ (ref. 30). To the order of accuracy of other assumptions made in this paper, it is reasonable to assume that z ∞ is sufficiently deep so that the temperature gradient can be taken as zero there; hence we take

As a robustness check, we computed results for small variations from ϕ=0, and found only minor quantitative differences between them, thus showing that this is a reasonable condition to impose.

To prescribe the underlying equations within each region, we need to consider the conservation of heat energy. While one could use the full heat equation (ρcT t =kT zz , where c is the specific heat capacity of the medium in question), we are able to simplify matters by considering the time scales involved. By noting that the critical depth scale in the model is the annual ice uplift height H (H=V × 1 yr), it allows us to compute the associated time scale of sinking Λ 1 , as ∼10 days (for the Antarctic parameters given in Table 2). The latter is found by comparing the heat flux required for sinking the meteorite with the solar forcing felt on it , and so . In contrast, the time scale Λ 2 for which the full heat equation will relax to its steady-state version (T zz ≈0), can be shown to be the order of 1 h (Λ 2 =ρ i cH2/k i ). With this large relative difference between time scales , we need only consider the steady-state version within each region j=1–4 (where we are implicitly assuming z ∞ is suitably shallow for the steady-state heat equation to be used, and simultaneously deep enough for the annual temperature variation to be minor), namely:

It is of note that time dynamics are still present within this model, via the annual variation in incoming solar radiation, thereby making our model quasi-steady. In addition, the steady-state heat equation coupled with the prescription of ϕ=0 for the Antarctic situation, sets the ice-temperature within region 4 as that of the meteorite base. This removes any dominating heat flux within region 4, making results computable from regions 1–3 only, and so an exact position for z ∞ is not required in that instance.

Rearrangement of the above equations yields a numerically tractable set of non-linear differential equations for a and b in terms of the input parameters. The particular form of their solutions depends on which stage of the sinking process is currently in effect. During the first stage, when the meteorite is encased in ice and moving with the upwelling ice, the speed of the interfaces a and b will be −V and the surface temperature of the meteorite will be below zero. In this case, calculation of the linear temperature profile is straightforward from equations (1, 2, 3, 4, 5, 6, 7), with no time evolution of a and b to solve for.

During the second stage, where melting of the lower ice has not yet commenced but the upper water layer is in existence (b−a>0), one can compute the location of the ice/water interface, a, from the dynamic equation

where

This last term in S net is necessary to ensure conservation of total solar radiation. Within this relation we are (in effect) assuming that any scattered shortwave energy within the ice only affects the ice at the atmospheric interface, that is, the shortwave energy directly warms only the ice surface and the meteorite (which then heats up the remaining ice by conduction). This assumption, which greatly simplifies the analysis, is expected to lower very slightly the temperature of the ice near the meteorite compared with reality. This consequently reduces the rate of sinking, and thus is, by design, a slightly conservative estimate (note that the numerical experiment confirms that, even when the last term in S net is removed completely, there is still little quantitative change to the results).

Once melting of the lower meteorite/ice surface (z=b+w) has commenced, and thus the meteorite starts to sink, one must switch to computing the coupled pair of equations,

To determine whether or not melting of the lower surface has begun, one can simply check whether or not the velocity of b in (14) is positive: while it is positive, we solve for a and b from (13) and (14), and if its numerical value ever becomes negative (no melting), one takes b as fixed and reverts to determining the evolution of a from equation (8).

The final aspect to note in the Antarctic representation of the model is that the atmospheric energy fluxes are seasonally dependent. To account for this we calculated and used the six-hourly average shortwave flux S(t) from the libRadtran atmospheric radiative-transfer model31, which incorporated a pseudo-spherical approximation with parameter inputs made relevant to the (relatively) well-parameterised Frontier Mountain Meteorite Trap area climatology32,33,34,35,36 (see Table 1 for details of meteorite collection in this locality). The period of this time average was chosen so as to maintain as much diurnal granularity as possible, without violating our assumption of a quasi-steady heat model, that is, 6 h Λ 2 ≈1 h (although results for durations between 1 h and 12 h showed only a minor quantitative difference). The computed maximum, minimum and mean daily values for the incoming shortwave solar flux S are shown in Fig. 4, where the diurnal variations are sinusoidal (and so the six-hourly averages lie within the shown range). To reflect the fact that the seasonality in shortwave energy has a direct effect on air temperature T a and thus the incoming longwave energy flux Q long , we allow these parameters to also be seasonally dependent37. To achieve this in a straightforward manner, we use a normalized version of the six-hourly shortwave profile, denoted by , to approximate their dynamics, namely:

Figure 4: Ice surface solar energy. The computed incoming shortwave energy flux S(t) reaching the Frontier Mountain meteorite trap ice surface, showing the daily mean (solid line), maximum daily value (upper dashed line) and minimum daily value (lower dashed line). These were calculated using the libRadtran atmospheric radiative-transfer model31 (as detailed in the Methods section). Full size image

where is a negative constant measured in °C (thus the maximum air temperature is attained during summer and the minimum temperature during winter), and and are both positive constants. With this model formulation, we are able to solve for the Antarctic situation, as well as the laboratory case (where θ=90° and the incoming energy fluxes are held constant).

The parameter values used in the laboratory simulation and the Antarctic simulation are all stated in Table 2. The in situ Antarctic parameter values are all in relation to the Frontier Mountain Meteorite Trap area19. Suitable proxies were taken for the blue-ice thermal conductivity k i (ref. 38), the surface roughness estimate u* ref. 30 and the solar attenuation parameter through water γ w (ref. 39) and ice γ I (refs 12, 39). When a range of parameters were provided, we used an averaged value. The percentage incoming solar radiation that is shaded, s h , is a time-averaged value inferred from the neighbouring mountain elevations17. The six-hourly solar elevation angles θ were computed40, while in the laboratory the solar-simulator lamp was held directly overhead (that is, 90°). The meteorite width w used in the Antarctic situation was chosen so as to be indicative of typical collected Antarctic meteorite specimens2 (although these smaller collected sizes/masses may be the consequence of larger material that broke up on the surface through repeated freeze–thaw cycles and wind action2). The scaling of the Antarctic longwave radiation was taken to be consistent with seasonal observations37. The thermal conductivity of the IIAB iron (Sikhote-Alin) was inferred from the IAB iron meteorite Campo del Cielo17; ±25% variations to this reasonable estimate still yielded results close to the data points (which was to be expected over the shorter laboratory time scales). The meteorite surface albedos (fusion crust cover) were independently measured for this study (see the Acknowledgements section).

The results were computed using a code written in Matlab, which the authors can supply on request. Potential extensions to our modelling approach are discussed in Supplementary Note 3.

Laboratory experiments

Our experiments centred around subjecting a meteorite encased in a block of ice (400 × 400 × 50 mm) to the radiation from an Oriel Solar Simulator arc lamp (irradiance 1,440 W m−2) that was held at a constant 10 mm above the ice surface and focused onto the englacial meteorite. This was conducted in a (otherwise dark, non-reflecting) temperature-controlled room, with ambient temperature −1 °C. The meteorite movement was recorded over a 3-h period using an HD time-lapse camera positioned at the side of the ice block, from which the meteorite’s progress could be measured (with an error of under ±1 mm), see Fig. 5, and to confirm that the outward-facing ice surfaces were not melting. With this experimental set-up, we were able to successfully conduct four controlled laboratory experiments, as shown in Fig. 1.

Figure 5: Laboratory stills. Laboratory measurements of an iron meteorite sinking through ice as discussed in the Methods section. Depth (z-axis) in mm is taken relative to the upper ice surface. Shown are three side-view images of the meteorite sinking downwards, taken at 80-min intervals. The images relate to the first data point (0 min), fifth data point (80 min) and ninth data point (160 min) of Fig. 1 (right panel), of the meteorite starting from 14.9 mm encapsulation depth in the ice. To align these images with the results of Fig. 1 (which show the depth of the upper meteorite surface) one must subtract the meteorite width of 10 mm from the depth of the meteorite base, which are indicated by the white lines (the location of the base is easier to observe, due to the reduced amount of glare from the light source). Full size image

To create the blocks of ice with a meteorite included, we first slowly froze the lower half the block in a container. We then placed a meteorite on top and carefully added an upper layer of cold water, which we then slowly froze. Once frozen, we removed the container. This whole process could take 2–3 days, so as to reduce the number of air bubbles within the ice. It also helped fully bond the meteorite to its surrounding ice.

Further to our four main experiments, controlled experiments were performed to show that the meteorite only sank in the presence of the solar simulator. We also recorded temperatures within the ice to confirm that when the meteorite was sinking its base was at 0 °C, thus confirming the assumptions within our energy balance model.