Biolocomotion and premelting in ice

Biota are found in glaciers, ice sheets and permafrost. Ice bound micro-organisms evolve in a complex mobile environment facilitated or hindered by a range of bulk and surface interactions. When a particle is embedded in a host solid near its bulk melting temperature, a melted film forms at the surface of the particle in a process known as interfacial premelting. Under a temperature gradient, the particle is driven by a thermomolecular pressure gradient toward regions of higher temperatures in a process called thermal regelation. When the host solid is ice and the particles are biota, thriving in their environment requires the development of strategies, such as producing exopolymeric substances (EPS) and antifreeze glycoproteins (AFP) that enhance the interfacial water. Therefore, thermal regelation is enhanced and modified by a process we term {\em bio-enhanced premelting}. Additionally, the motion of bioparticles is influenced by chemical gradients influenced by nutrients within the icy host body. We show how the overall trajectory of bioparticles is controlled by a competition between thermal regelation and directed biolocomotion. By re-casting this class of regelation phenomena in the stochastic framework of active Ornstein-Uhlenbeck dynamics, and using multiple scales analysis, we find that for an attractive (repulsive) nutrient source, that thermal regelation is enhanced (suppressed) by biolocomotion. This phenomena is important in astrobiology, the biosignatures of extremophiles and in terrestrial paleoclimatology.


I. INTRODUCTION
Ice sheets are an essential reservoir of information on past climate and they contain an important record of microorganisms on Earth, recording ice microbes and their viruses over long periods [1,2]. In these extreme environments, the abundance of virus is well correlated with bacterial abundance, but is 10 to 100 times lower than in temperate aquatic ecosystems [3]. Even in these harsh conditions, the virus infection rate is relatively high [4], leading to the expectation of low long-term survival rates. However, recent studies have shown that some viruses develop survival strategies to maintain a long-term relationship with their hosts [4,5], possibly up to thousands of years [6]. For example, viruses such as bacteriophages can switch to a lysogenic life strategy enabling them to replicate and maintain themselves in the bacterial population without lysis over multiple generations [4]. Moreover, among these viruses some can provide immunity to their hosts against other viruses [4,7], or manipulate their metabolism to facilitate nutrient acquisition by affecting motility genes [6]. Indeed, motile biota are found to be active in ice for substantial periods. For example, recently a 30,000 year old giant virus Pithovirus sibericum was found in permafrost [8] along with microbes [9,10] and nematodes [11], and viable bacteria have been found in 750,000 year old glacial ice [12]. Basal ice often contains subglacial debris and sediment, which serve as a source of nutrients and organic matter, providing a habitat for micro-organisms adapted to subfreezing conditions [13,14]. Additionally, the microbiomes of sediment rich basal ices are distinct from those found in glacial ice and are equivalent to those found in permafrost [13], expanding the nature of subfreezing habitats.
Ice cores provide the highest resolution records of past climate states [e.g., [15][16][17][18][19][20]. Of particular relevance to our study is their role as a refuge for micro-organisms, from the recent past [21,22] to millennia [23][24][25]. Ice microbes are taxonomically diverse and have a wide range of taxonomic relatives [14,24,[26][27][28]. Common algae taxa are centric and pennate diatoms, dinoflagellates and flagellates [29][30][31], whereas common bacterial taxa are pseudomonadota, actinobacteria, firmicutes and bacteroidetes [6,32]. Many of these microbes have different motility mechanisms [33,34] from swimming (e.g., Chlamydomonas nivalis [35] or Methylobacterium [6,36,37]) to gliding (e.g., diatoms [38,39] or Bacillus subtilis [24,40]), which can be used to assess their locomotion. Examples of biological proxies include diatoms [41] and bacteria colonies [42,43], reflecting a unique range of physical-biological interactions in the climate system. Therefore, understanding the relationship between the evolution of ice bound micro-organisms and proxy dating methods is a key aspect of understanding the covariation of life and climate. Finally, such understanding is essential for the study of extraterrestrial life. In our own solar system, despite the The nutrient source is shown by the purple gradient. The external temperature gradient induces a drift velocity, and particles move toward regions of higher temperature, in a process known as thermal regelation (black arrow). An additional drift velocity is associated with particle motion towards higher concentration of nutrients (black dotted arrow). Thus, depending on a particle's position, and the background temperature and nutrient gradients, these two drift effects can compete or amplify each other.
debate regarding the existence of bulk water on Mars [e.g., 44], thin water films, such as those studied here, hold the most potential for harboring life under extreme conditions. Indeed, lipids, nucleic acids, and amino acids influenced by active motility may serve as biosignatures of extra terrestial life. Combining measurements of diffusivity-characterizedmotility [45,46] with bioparticle distribution observed on Earth, provides crucial information for development of new instrumentation to detect the presence of extra terrestrial life [46,47]. Indeed, recently micro-organisms trapped in primary fluid inclusions in halite for millions of years have been discovered [48], providing promise for both terrestrial and extraterrestrial biosignature detection. When a particle is embedded in ice near the bulk melting temperature, the ice may melt at the particle-ice surface in a process known as interfacial premelting [49]. The thickness of the melt film depends on the temperature, impurities, material properties and geometry. A temperature gradient is accompanied by a thermomolecular pressure gradient that drives the interfacial liquid from high to low temperatures, and hence the particle migrates from low to high temperatures in a process called thermal regelation [49][50][51][52][53]. Thermal regelation of inert particles plays a major role in the redistribution of material inside of ice, which has important environmental and composite materials implications [49][50][51][52][53]. Moreover, surface properties are central to the fact that extremophile organisms in Earth's cryosphereglaciers, sea ice and permafrost-develop strategies to persist in challenging environments. Indeed, many biological organisms secrete exopolymeric substance (EPS) [54] or harness antifreeze glycoproteins (AFP) [55,56] to maintain interfacial liquidity. For example, sea ice houses an array of algae and bacteria, some of which produce EPS to protect them at low temperature and high salinity [57,58]. Additionally, the enhanced liquidity associated with high concentrations of EPS alters the physical properties of sea ice and thereby play a role in climate change [59,60].
In bulk solution, active particles act as simple microscopic models for living systems and are particularly accurate at mimicking the propulsion of bacteria or algae [e.g., [61][62][63][64][65]. By converting energy to motion using biological, chemical, or physical processes, they exhibit rich collective emergent motion from ostensibly simple rules [e.g., 66,67]. Algae and bacteria operate in complex geometries and translate environmental conditions into microscopic information that guides their behavior. Examples of such information include quorum sensing (e.g., particle population density), used by bacteria to regulate biofilm formation, defense against competitors and adapt to changing environments [68][69][70]; chemotaxis (e.g., concentration gradients of nutrients), used by algae/bacteria to direct their motion toward higher concentrations of beneficial, or lower concentrations of toxic, chemicals [71][72][73][74][75]. It is important to emphasize that factors such as surface adhesion, salinity, the segregation of impurities of all types from the ice lattice, among other factors [See e.g., 74,[76][77][78], make our treatment of chemotaxis a simplified starting point. However, field samples and laboratory experiments have shown that cell motility is influenced by chemotaxis at low temperature [45,79,80]. Thus, although there are many complicated interactions that provide scope for future work, the basic role of chemotaxis in the distribution of biota in ice must start with a self-consistent framework, which is the focus of our work.
The confluence of thermal regelation, bio-enhanced premelting and intrinsic mobility underlie our study. Indeed, intrinsic mobility and chemotaxis may compete with thermal regelation, which constitutes a new area of research-ice bound active particles in premelting ice, as illustrated in Fig. 1. Moreover, including micro-organism protection mechanisms that enhance interfacial liquidity, such as the secretion of EPS, constitute a unique class of regelation phenomena. Finally, treating this corpus of processes quantitatively is particularly relevant for climatology and the global carbon cycle [e.g., 81,82].
Our framework is the active Ornstein-Uhlenbeck particle (AOUP) [e.g [83][84][85][86][87][88][89]. The active force is governed by an Ornstein-Uhlenbeck process with magnitudeD a , which is the active diffusivity. This force can be compared to a colored noise process [85,90]. In addition to the active diffusivity, the AOUP is characterized by a time τ a , which defines the noise persistence, from which the system switches from a ballistic to a diffusive regime. The active diffusivityD a and characteristic time τ a can be measured experimentally [91][92][93]. The AOUP has been shown to provide accurate predictions for a range of complex phenomena [84,86,87,94], and is theoretically advantageous due to its Gaussian nature [85]. These issues motivate our use of the AOUP model framework to describe the motion of active particles in ice under an external temperature gradient with a nutrient source. We analyze these particles in three dimensions using a multiple scale expansion to derive the associated Fokker-Planck equation.
The paper is organized as follows. In order to make our treatment reasonably self-contained we note that we are generalizing our previous approach [53,95], which we recover in the appropriate limit. Thus, in §II we introduce the active Ornstein-Uhlenbeck model for bio-premelted particles and in §III we derive the associated Fokker-Planck equation using a multiple scale expansion. We then compare our analytic and numerical solutions after which, in §IV, we draw conclusions.

II. METHODS
Thermal regelation is understood as a consequence of the premelted film around a particle, originally treated as inert, that (a) executes diffusive motion in the ice column with diffusivityD(z)I, where I is the identity matrix, and (b) experiences a drift velocityṽ(z) = U (z)ẑ parallel to the temperature gradient [53]. Therefore, regelation biases the motion of an active particle by the drift velocity U (z)ẑ.
For inert particles with a sufficiently large number of moles of electrolyte impurities per unit area of surface, N i , the premelted film thickness d ∝ N i [53,96]. However, the production of EPS/AFP enhances liquidity at the ice surface by increasing the impurity concentration [14,59,77,97], which we treat here using an enhancement factor as N = nN i , which gives where the universal gas constant is R g , the latent heat of fusion per mole of the solid is q m , the molar density of the liquid is ρ l , the undercooling is ∆T = T m − T with T m = 273.15K the pure bulk melting temperature and T the temperature of ice. The velocity and premelting-controlled diffusivity are given by 6νRTm , with |∇T | the external temperature gradient. The viscosity of the fluid is ν, the particle radius is R and k B is the Boltzmann constant. Here, ρ s q m ∼ 334 × 10 6 J m −3 [53]. The evolution of the particle positionr = (r 1 ,r 2 ,r 3 ) = (x,ỹ,z) and its activity are described by two overdamped Langevin equations The first term on the right-hand side of Eq. (4) treats the chemotaxis response, representing the effect of the nutrient source of concentrationC on the particle dynamics, where β D is the chemotactic strength [98][99][100][101], which we treat as a constant determined by the parameters in our system. We note, however, that the transport properties of sea-and glacial-ice depend on their unique phase fraction evolution [e.g. [102][103][104][105], which would clearly influence the effective-porosity dependent β D . In the ideal case, wherein the nutrient source is isotropic and purely diffusive, whereD ch is the nutrient diffusivity. The activity, or self-propulsion, is given by the term 2D aη in Eq. (4), with D a the active diffusivity. The latter represents the active fluctuations of the system, such as those originating in particular processes described in Refs. [106][107][108][109]. Nutrient sources, such as dissolved silica, oxygen, nitrogen and methane, play a vital role in the life of ice-bound micro-organisms, such as algae and bacteria [e.g., [110][111][112][113][114][115]. Here we assume thatD ch >D a , consistent with [116][117][118], andD ch >D(z). The functionη = (η 1 ,η 2 ,η 3 ) is described by an Ornstein-Uhlenbeck process, with correlations given by where τ a is the noise persistence as noted above. In the small τ a limit,η reduces to Gaussian white noise with In contrast,η does not reduce to Gaussian white noise when τ a is finite, and Eq. (4) does not reach equilibrium. Hence, τ a controls the non-equilibrium properties of the dynamics [85,86]. Finally, the random fluctuations in Eqs. (4) and (5) are given by zero mean Gaussian white noise processes The Langevin Eqs. (4) and (5), allow us to express the probability of finding a particle at the positionr = (r 1 ,r 2 ,r 3 ) = (x,ỹ,z) at a given timet through the Fokker-Planck equation, which describes the evolution of the probability density function P (r,η,t|r 0 ,η 0 ,t 0 ), with the initial condition P (r,η,t =t 0 |r 0 ,η 0 ,t 0 ) = δ(r−r 0 )δ(η−η 0 ). To simplify the notation, we write the conditional probability as P (r,η,t|r 0 ,η 0 ,t 0 ) ≡ P (r,η,t) and eventually arrive at the following system of coupled equations, Equations (8) and (9) describe the space-time evolution of the probability of finding a particle and the concentration of nutrients respectively, akin to those of [83,85,119], but including the effects of thermal regelation discussed above. Both equations contain microscopic and macroscopic scales. The regime of interest is the long time behavior, computed by deriving the effective macroscopic dynamics as described next.

A. Method of multiple scales
The macroscopic length characterizing the heat flux is The particle scale l is such that l << L, and hence their ratio defines a small parameter = l L .
We use the microscopic length l and a characteristic time τ , determined a posteriori, and introduce the following dimensionless variables whereṽ a = 2Da τa [120,121], v ac is the characteristic active velocity, u and D c are the characteristic values of the regelation velocity and the premelting enhanced diffusivity respectively, and D n and c h are the characteristic values of the diffusivity and nutrient concentration respectively. With these scalings, Eqs. (8) and (9), become and (14) in which we have the following dimensionless numbers, We identify four characteristic time scales: t diff l = l 2 /D c , t adv l = l/u, t diff L = L 2 /D c and t adv L = L/u, associated with "microscopic" diffusion and advection on the particle scale, l, and "macroscopic" diffusion and advection over the thermal length scale, L. Nutrient and premelting enhanced diffusivity are taken to operate on the same time scale; t diffn l,L ∼ t diff l,L . The Péclet number represents the ratio of the characteristic time for diffusion to that of advection, and those associated with regelation and activity are P e and P a respectively, and can be defined over both length scales, The temperature gradient across the entire system drives thermal regelation and hence advection dominates on the macroscopic scale, so that P L e = O(1/ ), or equivalently, t adv . On the macroscopic scale P L e becomes large, as tends to zero, and thus we use the macroscopic advection time τ = t adv L as our characteristic time, so that Péclet numbers based on L are O(1/ ) and those based on l are O( ). In consequence, Eqs. (13) and (15) Now, we let R =r/L describe the macroscopic length scale, and T =t/t adv l describe the microscopic time scale, leading to the following stretching of the microscopic scales; Next, we use a power series ansatz for the state variables, to derive a system of equations at each order in [122], which for the concentration of nutrients, Eq. (19), are shown to second order. We take the approach described in [123,124] to solve Eqs. (23)- (25). We integrate by parts over the microscale variables r and use the periodic boundary conditions to obtain the so-called weak formulation [125] of the leading order Eq. (23), the solution of which relies on the following product ansatz C 0 (r, R, T, t) = ζ(r)c 0 (R, T, t) . The existence and uniqueness of C 0 is ensured using the Lax-Milgram theorem [125], also known as the solvability condition or the Fredholm alternative [123]. Thus, C 0 is constant over C 0 (r, R, T, t) = C 0 (R, T, t). The solvability condition for the equation at O( ) is from which we find that c 0 is stationary over T , leading to C 0 (R, T, t) = C 0 (R, t) and C 1 (R, T, t) = C 1 (R, t). Substituting C 1 into the O( 2 ) Eq. (25) and using the solvability condition, gives nutrient diffusion on the macroscale as showing that, as expected, the homogenization procedure is consistent with the well-known self-similar behavior of diffusion [126]. The order by order equations for the probability density function described by Eq. (18) are simplified by the observation that C 0 and C 1 do not depend on r, and C 0 only contributes at order O( 2 ), and hence we obtain where L = M + Q, with M = − ∂ ∂r3 v − P a v a η · ∇ r + ∇ 2 r D, and Q = P A ∇ η · η + P A 2 ∇ 2 η . Finally, as shown in Appendix V A, upon substitution of P 1 into Eq. (31) and using the solvability condition, we obtain the effective macroscopic dynamics as which are the dimensional forms of Eqs. (31) and (28) respectively. These capture the long time behavior wherein the active force is treated through the effective diffusivity, which is enhanced by thermal regelation, consistent with our previous work [95] and that in active matter systems generally [62,84,91]. Eqs. (32)- (33) can be mapped onto the well-known Keller-Segel equations for chemotaxis [99][100][101]127], where ρ is the cell density and the sign of β D determines whether a cell is attracted or repelled by the nutrient. Finally, when nutrients are neglected, β D = 0, we recover our previous results [53,95]. Although Eq. (32) has an analytical solution in the large Péclet number limit, which previously allowed us to study the effect of the activity ([e.g., see 53,95] or Appendix V B), here we fix the activity and focus on the competition between thermal regelation and bio-locomotion that require solving Eqs. (32) and (33) numerically. We show dimensional results because of our specific interest in these processes in ice.
In order to study the effect of nutrients on bio-locomotion, we fix the interfacial concentration of impurities and vary the chemotaxis strength β D , where nutrients either attract (β D > 0) or repel (β D < 0) the bio-particles. Because we are interested in the situation wherein the effects of chemotaxis compete with regelation, this constrains the magnitude of β D as follows. We ask for what order of magnitude of β D are the typical chemotactic speeds approximately the same as the regelation velocity in Eq. (4). Figure 3 shows the Gaussian solution of the concentration field, with a flux that becomes arbitrarily small at long times, dominated by the algebraic contribution to ∇rC(r,t) ∝ xt −3/2 exp(− x 2 t ) ∼ xt −3/2 . For the parameters studied here, the regelation speeds are 10 −12 − 10 −10 m s −1 [53,95], and hence we capture this same range in β D ∇rC(r,t), with |β D | = 10 −10 m 2 M −1 s −1 , which is realized across a large time span wherein ∇rC(r,t) varies by several orders of magnitude. This is also reflected in the dimensionless ratio β D c h Dc in Eq. (13). Namely, for micron to nm scale premelted films surrounding micron sized particles D c ranges from about 10 −14 −10 −13 m 2 s −1 , and the nutrient concentration over relevant time scales has mean values ranging over 10 −3 −10 −2 M. Therefore, β D c h Dc ranges from 1 to 100 and hence chemotaxis is on a similar footing to regelation under these circumstances. For all cases considered here we use β D = ±10 −10 m 2 M −1 s −1 for attractive/repulsive chemotaxis.
For attractive chemotaxis (β D > 0), we show in Figure 4 (a) the dependence of ρ(r,t) along thez-axis parallel to the temperature gradient and atx =ỹ = 0, with the concentration of nutrients centered atz = 55m. For the same conditions in the absence of chemotaxis, the net displacement from low to high temperatures due to regelation is approximately 10 m [95]. We see here the chemo-attractive modulation of ρ(r,t) during this displacement, which "pulls up" the high temperature (lowz) tail towards the lower temperature (largez) but higher concentration regions centered atz = 55m. The associated asymmetry depletes/attracts the low temperature regions at largerz and concentrates the high temperature regions at smallerz, and is reflected in the evolution towards a sigmoidal region transecting the source atz = 55m. As the maximum of ρ(r,t) advects through the source region it first sharpens, due to the chemo-attraction from the source "behind" it atz = 55m, and then begins to spread out again because of the decay in the chemotactic gradient in time as seen in Figure 3 and discussed above.
For repulsive chemotaxis (β D < 0), we see in Figure 4 (b) the broad sharpening of the initial distribution in the lower temperature (largez) regions as it regelates/advects into the diffuse repulsive tail of nutrient field to the right of the source region centered atz = 55m. However, because the initial high temperature (smallz) tail of ρ(r,t) interacts with the nutrient source region atz = 55m, chemo-repulsion quickly drives particles towards high temperature (smallz) regions, and is clearly reflected in the creation of a local maximum. This maximum advects towards high temperature with a decaying amplitude and width due to the decay in the chemotactic gradient in time as seen in Figure 3.
In Fig. 5, we show the combined effects of EPS/AFP surface enhancement of impurities in the absence of chemotaxis (β D = 0), as shown in Fig. 2, and the influence of chemotaxis on particle dynamics for fixed surface impurities, as shown in Fig. 4. As we vary the surface concentration of impurities we observe the same basic features as described in Figs. 2 and 4 and hence the same physical description applies in their interpretation. Namely, regardless of whether chemo-attraction or chemo-repulsion is operative, if the interfacial concentration of impurities N is sufficiently large then the interfacial film thicknesses are sufficiently thick that thermal regelation dominates the evolution of ρ(r,t).
As the interfacial concentration of impurities N decreases chemotaxis exerts more control on the distribution, and the basic dynamics are the same as described in Fig. 4. Because the magnitude of β D is fixed, and the characteristic concentration c h is 10 −2 M, this N -dependence is simply assessed as discussed above, in terms of the dimensionless ratio β D c h Dc in Eq. (13). Namely, the numerator is fixed, but as N increases so too is the film thickness d through Eq. (1), and since D c ∝ d 3 [52], then β D c h Dc decreases as N −3 , and the balance between chemotaxis control of the distribution gives way to regelation control. The corpus of effects studied here are reflected in this basic balance and shown in Figs. 2, 4 and 5.

IV. CONCLUSION
Micro-organisms in ice exhibit complex processes to persist and evolve in their harsh environments. They have developed different survival strategies, such as producing exopolymeric substances or antifreeze glycoproteins, and directing their motion toward nutrients or away from waste [39,74,132,133]. We have modeled such micro-organisms using active Ornstein-Uhlenbeck particles subject to thermal regelation and biolocomotion in three dimensions. Firstly, we used a multi-scale expansion to derive the relevant coupled Fokker-Planck and diffusion equations (32)- (33). Secondly, when nutrients are neglected, and the chemotactic strength β D = 0, we model the bio-production of surface chemicals, such as exopolymeric substances or antifreeze glycoproteins, as a surface colligative effect, and find that the associated bio-enhanced thermal regelation can dominate the distribution of particles in ice. Consistent with previous results [95], in a large Péclet number limit analytical solutions for the particle distributions are possible, and are consistent with the numerical solutions as shown in Fig. 2. Thirdly, we studied the competition between thermal regelation and biolocomotion, as function of the chemotaxis strength β D , the interplay between which is shown in Figs. 4-5. The relative importance of chemo-attraction and chemo-repulsion to thermal regelation is captured by the dimensionless ratio β D c h Dc . When this ratio is large we find a complex modulation of regelation by chemotaxis, and when small, due to increased surface impurity concentration, leads to regelation dominated redistribution of particles. We note, however, that we have not treated the process wherein nutrients themselves have a colligative effect, which would introduce a particularly complex spatio-temporal dynamics.
Finally, we describe settings to which our analysis is applicable. It is of general interest to understand how particles in ice migrate in response to environmental forcing, as they are used as proxy to infer past climate [19,134,135]. Moreover, bioparticles in ice migrate in response to environmental forcing, and micro-organisms play an important role in climate change [136][137][138]. For example, an increase in temperature activates algae/bacteria trapped in ice, producing chemicals that increase their mobility [138]. Indeed, an increase in algae/bacteria decreases the albedo of the ice [139][140][141], thereby enhancing melting. Finally, understanding the distribution and viability of bioparticles in partially frozen media on Earth [e.g., 142,143] is essential in astrobiology [47,144,145].

CONFLICT OF INTEREST STATEMENT
The authors declare no competing interests.   The solution of the leading order Eq. (29) is derived by making the following product ansatz P 0 (r, R, η, T, t) = w(r, η)ρ 0 (R, T, t) .
Integrating by parts over the microscale variables r and η, and using the solvability condition, the solution of Eq. (29) P 0 is constant over the period: P 0 (r, R, η, T, t) = P 0 (R, η, T, t). The leading order Eq. (29) becomes Using the known result for a multi-dimensional Ornstein-Uhlenbeck process [146], the solution for w is given by The solvability condition for O( ) equation is which depends on the leading order result, P 0 , from which we find and the O( ) equation becomes We assume that after which we find that Substitution of P 1 into the O( 2 ) equation and using the solvability condition, we obtain and in dimensional form, we have B. Solution of the Fokker-Planck equation in the large Péclet limit The analytic solution of Eq. (32) in the large Péclet number limit follows by expanding the probability density function, ρ perturbatively [147], as with ρ 1 = O(c)ρ 0 . The associated system of equations is To derive the solutions for ρ 0 and ρ 1 , we use the Green's function G(r,t|r 0 ,t 0 ) [148] which satisfies The solution of Eq. (47) was derived in the large Péclet number limit in [95], and is ρ0 = G(r,t|r0,t0) =z .