Near-perihelion activity of comet 67P/Churyumov–Gerasimenko. A first attempt of non-static analysis

The observed rate of water production of comet 67P/Churyumov–Gerasimenko near its perihelion can be approximated by a very steep power function of the heliocentric distance. Widely used thermophysical models based on a static dust layer on top of the icy/refractory matrix are poorly consistent with these observations. We analyse published model results and demonstrate that thermophysical models with a uniform and static ice free layer do not reproduce the observed steep water production rates of 67P near perihelion. Based on transient thermal modeling we conclude that the accelerated gas activity can be explained assuming that the active area fraction near the south pole is increased. The deeper penetration of the heat wave during polar day (no sunset) can activate sublimation through thicker inert dust layers. This can also lead to removal of thicker dust layers and consequently to an expansion of the


I N T RO D U C T I O N
Comet 67P/Churyumov-Gerasimenko surprised the Rosetta scientists with its complex shaped nucleus. It is covered with steep cliffs, some several hundred metres high, with cavities and pits, with rugged terrain predominantly on the Southern hemisphere and with flat smooth dust covered surfaces (Thomas et al. 2015). The bi-lobate nucleus rotates with an obliquity of 52 • (RA = 69 • and Dec. = 64 • ), such that the southern solstice coincides with its perihelion passage, and a pre-perihelion period of 12 h 40 m (Mottola et al. 2014). The size of the nucleus is about 2 × 4 × 3 km with a volume of 18.56 ± 0.02 km 3 (Preusker et al. 2017). Southern summer coincides with the perihelion passage and lasts only 8 months out of the 6.4 yr orbital period while the comet travels from the pre-perihelion equinox at R h = 1.6 to R h = 1.8 au postperihelion. Typically erosion due to sublimation of water ice is four times stronger on the Southern hemisphere than on the northern areas. Only a fraction of the surface has to be active to match the observed peak activity during perihelion (Keller et al. 2015b). The activity of the nucleus (dust and water production) appears to be rather uniform, controlled by insolation in spite of the widely varying heterogeneous morphology. A simple two layer (a thin inert dust layer on top of the cometary dust/ice matrix) thermophysical E-mail: skorov@mps.mpg.de model, assuming homogeneous activity, was successfully applied to calculate the torque that changes the rotation period of 67P (Keller et al. 2015a).
This simple model, however, cannot reproduce the steep increase of the water production rate (absolute value of exponent >5) observed a few months before perihelion [see discussion in Attree et al. 2019) and models presented in Hansen et al. (2016) and Fougere et al. (2016)]. Forward modelling of the coma using harmonic expansion and DSMC method required spatial variability of activity to match the in situ density (ROSINA) and remote sensing observations (MIRO and VIRTIS) (see Combi et al. 2020) for an overview). Selected distributions of surface areas with different intensity of activity do influence the exponent of the production rate variations taking advantage of the complexity of the nucleus (Marshall et al. 2017). The same approach is required to simulate the observed non-gravitational forces . The approximations obtained are based on simple parametric fitting without analysis of the required microphysical constraints. All these models apply a static thermophysical approach where the surface properties (thickness of the dust layers) do not change with heliocentric distance. However, the seasonal variation of the surface colour from reddish to bluish when the comet approached the sun and back again to reddish when it receded (Fornasier et al. 2016) suggests that the surface properties change during the perihelion passage, e.g. that the surface becomes more ice rich (bluish) when activity level and insolation increase. Attree et al. Table 1. The exponent for the power functions approximating the total gas production (models and observations). A more detailed description of the models and corresponding graphics is presented in the Appendix. (2019) achieved the best fit to match the non-gravitational forces when they assumed that the surface activity strongly increased near perihelion, stronger than the insolation (see fig. 5 in Mottola et al. 2020). Around perihelion mainly the areas near the south pole are illuminated while they are unilluminated during northern summer. In this paper, we demonstrate that the deeper penetration of the heat wave during polar day (no sun set) increases the active area fraction near the south pole causing both the awakening of previously inactive sub-domains, and possibly the destruction of a weakening dust crust. These effects can be responsible for the observed excess productivity.

RO L E O F N U C L E U S S H A P E
We begin our analysis with a brief consideration of the simplest hypothesis: the accelerated growth of gas production a few months before perihelion can be caused not only by an increase in solar radiation absorbed by the unit area of the cometary nucleus but also by an increase in total illuminated area. This hypothesis is simple and easily verifiable. One of the options for such verification is presented in Marshall et al. (2019). Using the simplest model where gas production is calculated from the energy balance equation on the surface of a cometary nucleus (Model A in Keller et al. 2015b) and analysing gas activity as a function of heliocentric distance, the authors investigated the role of the nucleus shape and its spin axis orientation for several short-periodic comets. They found that in case of 67P, the illuminated cross-section does not significantly change over its orbit. As a consequence, the gas production of 67P varies as function of the heliocentric distance similarly to a nucleus of an appropriate spherical shape. Below we present the results obtained for Model A based on the most recent shape model of the nucleus (SHAP7; Preusker et al. 2017) decimated to 125 000 facets. We also add into the model shadowing and self-illumination effects.
As it was done in Keller et al. (2015b), for each orbital point the gas production is calculated explicitly taking into account the rotation of the nucleus at 50 azimuth positions. The obtained numerical results were approximated by a power function (R H < 3au). The exponent of this approximation is presented in Table 1. Besides we show results obtained for (1) a spherical nucleus with the same surface area, (2) a nucleus with the observed shape, but with an axis of rotation perpendicular to the orbital plane (zero obliquity) with self-illumination and without self-illumination. Accounting for self-heating and shadowing makes only minor changes to the water production of 67P inside ∼3 au. Note that in all tested cases, gas production does not show accelerated growth as observed after the comet passes through the inbound equinox (10 May 2015, 95 d before perihelion at a heliocentric distance of 1.67 au). All calculated curves are well approximated by a power dependence with a constant exponent (the correlation coefficient is very close to unity). 1 The obtained absolute values of the exponent are less than 3 varying between 2.3 and 2.6 that is significantly smaller than the observed values (Fougere et al. 2016;Hansen et al. 2016;Marshall et al. 2017;Läuter et al. 2019). motion. Therefore, it can be expected that the observed acceleration in gas production will be explained by these models.
The recent theoretical models mentioned above, all study comet 67P, but differ markedly in detail. The authors focus on various specific features of the transfer processes (e.g. weakening of the efficiency of sublimation due to dust layer resistance, radiative thermal conductivity that plays an important role at high temperatures, surface fraction of ice and bulk ice/dust ratio that change cooling effectiveness of sublimation). They use either an assumption about a thermal relaxation of the surface layer (Models B and C in Keller et al. 2015b;Keller et al. 2017;Skorov et al. 2017) and retrieve an effective gas production from the system of algebraic equations describing the energy balance at the surface and at the ice/dust interface, or solve a classical differential transient heat transfer equation in order to evaluate the temperature of sublimating ice and, hence, an effective gas production (Shi et al. 2016;Skorov et al. 2016;Blum et al. 2017;Hu et al. 2017) In Table 1, we present the results of the power approximation obtained for these two-layer models. Approximation curves and additional fitting parameters, as well as some details of the models can be found in the Appendix. An accurate comparison of the exponents of the power functions obtained for the various numerical models is rather arbitrary. Nevertheless, some qualitative observations can be made based on such an analysis. The layers considered in Keller et al. (2015b) are very thin (<1 mm) and rather opaque for the gas flow (dimensionless thickness is 10 and 20). Whereas the layers used in Hu et al. (2017) and Blum et al. (2017) are much thicker (from 5 to 10 mm). In Model B and the first variant from (Hu et al. 2017), the dimensionless thicknesses are the same (that is, the gas-dynamic resistances are the same) and the effective thermal conjunctivitis are small and, apparently, comparable. We see that the exponents are very close and remain less than 3. The thicknesses of these layers are comparable to the expected penetration depth of the daily heat wave. As the model layer thickness increases (see Table 1), the exponent increases noticeably, i.e. the slope of the gas production curve is getting steeper. This observation is important and will be needed in the future. It is more difficult to analyse the results obtained in (Blum et al. 2017). The authors considered a monolayer of ice-free dust (thickness equals a grain size). Such a layer is very transparent for a gas flow (see e.g. Gundlach, Skorov & Blum 2011). In addition, the size of the model particles is very large: ten times larger than the size used in Hu et al. (2017) and a hundred times larger than in Keller et al. (2015b). Such a size leads to the fact that the non-linear contribution of radiation thermal conductivity becomes dominant. Perhaps this explains the lower value of the exponent for the same layer thicknesses (10 mm). In this model, the significant changes in the ice fraction do not change visibly the slope of the gas production curve (see Fig. A4).Thus, we can conclude that in all the cases considered, the gas production after passing through the inbound equinox is well described by a power function with a constant exponent.
It is important to note that in all these models it is assumed that the surface properties (first of all, the thickness of the porous dust layer) remain unchanged both with a change in the heliocentric distance and for different regions of the nuclear surface. This assumption can be called a homogeneous static dust layer. Here, it is important to explain one non-trivial feature of this assumption: the concept of static does not mean that the layer consists of the same particles, quite the opposite, it is assumed implicitly that the dust particles on the upper boundary of the layer are constantly carried away by the gas flow, and this ablation is exactly offset by the addition of new particles on the lower boundary of the layer (so we have at the same time a fixed layer thickness and non-zero dust production). It seems obvious that this assumption is artificial and non-physical. However, today we do not have more adequate models. It is not surprising, therefore, that none of the models mentioned can fit an accelerated gas production near perihelion.

VA R I A B L E AC T I V I T Y
At the same time, it was shown (Keller et al. 2015b;Attree et al. 2019) that the observed effects are well simulated by the hypothesis of an ad hoc variable areal distribution of the ice fraction. As we noted earlier (see fig. 5 in Mottola et al. 2020), although the model presented by Attree et al. (2019) takes into account the different activities in the 26 named geographical regions of Thomas et al. (2015), their results also break down around perihelion, suggesting that temporal variations of activity are required in addition to the geographical ones. Therefore they invoked time-dependent activity in the Southern hemisphere, since the discrepancies between the observations and time-independent models were manifested around perihelion when the South dominates activity. Sharply increasing the Southern active fraction (by ∼3 times) in about three months before and three months after perihelion, they achieved good fitting of the non-gravitational forces. Unfortunately, in the proposed form, their hypothesis is essentially just a model pasteurization. Nevertheless, it gives us important clues about what processes should be included in future models. The question arises: How can we significantly increase activity in a microphysical model built on first principles?

P O L A R DAY A N D WA K E U P O F S L E E P I N G A R E A
To assess gas production, it is important to consider (1) how the energy required for ice sublimation is delivered to an icy region, (2) how the product of sublimation (vapour) is evacuated through a porous layer. The first process is controlled by the effective thermal conductivity, and the second by the permeability of the layer. In general, both of these characteristics are functions of the microscopic structure of the layer. In the simplest case discussed below, when a layer consists of quasi-spherical aggregates of one size, the transport characteristics of the layer are functions of the size of the aggregates R A and the ratio of its thickness L to this size (i.e. dimensionless layer thickness). The observations by Rosetta instruments do not allow us to significantly narrow the range of values of the thermal inertia of the surface layer. According to published results, its values cover an order of magnitude: from 10 to 100 JK −1 m −2 s −1/2 (Choukroun et al. 2015;Schloerb et al. 2015), i.e. the thermal conductivity is uncertain by a factor of hundred. The structural characteristics of the surface layer are even less known. This forces us to consider at the first stage these characteristics as free model parameters.
Hereinafter, we use a non-stationary thermal model in which the boundary condition for the ice/dust interface is set explicitly. A full description of the model can be found in Skorov et al. (2016). The main difference from the model used in Hu et al. (2017) and Shi et al. (2016) is that we take into account conductivity by thermal radiation, which plays a dominant role for relatively high temperatures and large particles (Gundlach & Blum 2012). The main idea related to the conductivity by thermal radiation is very simple: a sharp increase in thermal conductivity with increasing temperature (k R ∼ T 3 ) leads to an increase in the energy that is delivered deep into the nucleus and becomes available to ice Figure 1. Gas production as a function of heliocentric distance estimated for a spherical nucleus with orbital parameters, obliquity and rotation as for comet 67P. Results for high (solid curves) and low (dashed curve) thermal conductivity are shown. High thermal conductivity includes both solid and radiation thermal conductivities, whereas only the first type is taken into account in the low thermal conductivity model. Thick (50 mm) and thin (5 mm) layers are examined. Dimensionless thickness is 50. Positions of maximum production are marked. The time-step of the calculation is about 5 min.
sublimation. It can then be expected that corresponding accelerated growth in gas production will be observed at small heliocentric distances. Note that for calculating the resistance of a porous layer for gas flow in the Knudsen regime only the ratio of pore size (or particle size) to thickness is of importance. However, dimensionless scaling is unfortunately not applicable for calculating the thermal resistance of a dust layer. Therefore, below we have to consider layers having different thicknesses and different grain sizes. Note that based on the results presented above, we consider the spherical nucleus having an orbit, obliquity and rotation as 67P. Because we focus on changing overall gas production over heliocentric distance, this simplification is reasonable. Fig. 1 shows the gas production as a function of heliocentric distance obtained for two model layers (5 and 50 mm thick) having the same dimensionless thickness (thickness/grain radius) equal to 50. For a thin layer, variants with and without radiative thermal conductivity are shown. For a thick layer, the results are shown only for the case when radiation thermal conductivity is taken into account. The minimum value of solid thermal conductivity is about 0.003 W m −1 K −1 , and the maximum value of the total effective thermal conductivity is about 0.1 W m −1 K −1 .
As expected, in the case of low thermal conductivity, gas production near perihelion is almost an order of magnitude lower than gas production in a model with a high (solid + radiative) conductivity. The more unexpected result is obtained by comparing the results for thin and thick layers. Gas production from under a thick layer grows visibly faster: sublimation from under a fivecentimetre layer is more than ten times less than sublimation from under a five-millimetre layer 150 d before perihelion (i.e. near inbound equinox), and almost the same at perihelion. This effect is associated with a change in illumination behaviour. For a daily heat wave, the layer thickness is too large and only a small fraction of the absorbed solar energy is available for sublimation. For the southern subpolar region, the sun does not descend beyond the horizon (polar day). As a result, the heatwave penetrates deeper regions and the amount of energy available for sublimation increases rapidly.
If we now imagine that the surface is covered with layers of different thicknesses (in our model 5 and 50 mm), then as the comet approaches perihelion, the sleeping areas covered by the thick dust layer wake up. Their contribution to the total gas production increases and produces accelerated growth in qualitative agreement with the observations. Note that to produce this effect we need the presence of both thin and thick layers. If we exclude thin layers from the model, then the case of only a thick layer cannot match the observed activity at large heliocentric distances. We emphasize again: in our model, areas covered with a thick layer become active only because of the change in the illumination (the polar day effect).

P O L A R DAY A N D D E M O L I T I O N O F C RU S T
During polar day illumination, the heat wave penetrates thicker inert layers to enhance the gas production. Now we take the next step and investigate the question of whether an 'awakening' of a thick layer can lead to its demolition. If the thick layer is removed (and it has to be removed otherwise activity would be quenched) the area covered by a thin layer is increased and, hence, a gas production grows. Obviously, the condition for layer destruction is reduced to analysing the balance of gravity of the nucleus, cohesion of the material and the lifting force of the outgoing gas (see e.g. Gundlach et al. 2015;Skorov et al. 2017). The first force is constant in the approximation of a spherical nucleus. The second depends on the microstructure of the layer (Skorov & Blum 2012). The third force depends both on the microstructure of the layer (primarily on permeability) and on the energy available for sublimation. We will check whether (1) the pressure at the lower boundary of the thick layer is higher than the corresponding pressure for the thin layer, (2) this pressure is comparable or exceeds the force of cohesion. To make the model as clear as possible, we will consider layers with identical particles. The illustrating results are presented in Fig. 2.
We show the vapour pressure for two latitudes (equatorial and at 60 • of south latitude) and two layer thicknesses (5 and 50 mm). The aggregate size is 1 mm. Variants with and without the radiative conductivity are considered. Comparing the simulation results, we find that the maximum temperature under the thick layer is about 7 K higher than the corresponding temperature under the thin layer. This seemingly small difference in temperature leads to the fact that the corresponding pressures differ almost three times due to the exponential dependence of the saturated vapour pressure on temperature. We emphasize that the maximum pressure value under the thick layer (∼1 Pa) is higher than the theoretical estimate of the tensile strength of the layer under consideration (Skorov & Blum 2012). Thus, for the first time we show that for a reasonable set of model parameters, demolition of a dust crust due to the sublimation of water ice can happen on the surface of comet 67P. This result does not cancel our previous conclusion (Skorov et al. 2017), but modifies it: high thermal conductivity and polar day conditions, allowing the thermal wave to penetrate into deeper layers, are necessary for the crust removal. It is worth noting that the scenario is realized when the pressure under a thick layer overtakes and exceeds the pressure under a thin layer. This means that there is a time interval when a layer can grow remaining stable, and only after reaching a critical thickness, it is destroyed. The time span for the pressure under a thick layer to exceed the pressure under a thin one varies for different latitudes. Without making assumptions about the specific strength of the layer, this means that the destruction of a thick layers does not occur simultaneously but in a certain time interval. Note that a layer demolition always leads to an additional increase in gas production. Both, theoretical estimates (Skorov et al. 2017) and experimental results (Bischoff et al. 2019) suggest that the full pieces of a porous dust layer is removed in whole pieces, that is, there is a transformation from sublimation under the covering dust layer to the mode described by Model A: sublimation from an icy surface. As noted in Keller et al. (2015b), this model corresponds to the maximum gas production for any given illumination condition. The described qualitative scenario for the first time shows how a change in the structural properties of the covering dust layer could occur, and how these changes can lead to a significant acceleration of gas production near perihelion of 67P.

C O N C L U S I O N
The focus of our attention was the sharp increase in gas production observed after comet 67P passed the inbound equinox. The production rate increases much faster than R −2 H , with an exponent between −5 and −7 (Fougere et al. 2016;Hansen et al. 2016;Marshall et al. 2017;Läuter et al. 2019). As was already previously shown , this observation cannot be satisfactorily explained in the framework of the simplest model in which sublimation is calculated on the basis of the energy balance on a uniform surface.
Extending this analysis, we investigated the role of the nucleus shape, obliquity, and its rotation as well as shadowing and selfillumination effects. We have shown that in the assumption under consideration (the so-called Model A), gas production is very well approximated by a simple power function with an exponent less than three in all the tested cases. This result is qualitatively inconsistent with the observations noted above.
In the next step, we examined a class of more complex models. These models take into account the presence of an inert porous dust layer on the surface and the heat flux into the depth due to thermal conductivity. For these models, we investigated the influence of (1) gas resistance and thermal insulation of the icefree porous dust layer; (2) radiative heat conductivity that can be many times greater than the contact thermal conductivity for large dust particles; (3) a factor of 'icy area fraction', which also modifies the energy balance at the icy interface; (4) an inert layer of variable thickness. A factor of 'icy area fraction' is considered in detail in Hu et al. (2017) as a measure of the microscopic areal abundance of water ice at depth, where the dust and ice are in thermal equilibrium.
Summarizing the results, we conclude that the steep increase of the water production rate near perihelion of 67P cannot be matched by a two-layer model with a static porous dust layer covering the dust/ice matrix. At the same time, as it was shown by Attree et al. (2019) the observed effect is well matched by the hypothesis of a transient activity, i.e. a time-dependent active fraction, with the active fractions of two regions in the Southern hemisphere allowed to vary between two values. The ad hoc time dependent activity variation of the south pole region required to fit the non-gravitational forces ) can be explained by assuming the inert dust layer to be of variable thickness. Thicker layers are 'awakened' when the sun does not set anymore while the comet approaches southern solstice near its perihelion. Due to the permanent insolation, the heat wave penetrates to deeper levels. At the same time the pressure underneath thicker layers reaches higher values and can finally reach levels that can blow off the cover. The uncovered areas enlarge the region of normal activity. Keller et al. (2015b) showed that only 1/6 of the surface needs to be active needs to fit the observations (their Model B calculations). The observed roughness and morphology of the nucleus of comet 67P persuasively suggest that the thickness of the covering inert dust layer is variable due to varying insolation. Activating originally inactive areas can hence lead to a dramatic increase in gas and dust production, far beyond what is expected from the increase of insolation. The expansion of nucleus surface exposed to polar day conditions can easily explain the steep power exponents when the comet approaches perihelion (southern solstice).
In a further step, we will calculate the overall gas production rate by expanding the thermophysical model for a dust cover of variable thickness.

A P P E N D I X A : D I S C U S S E D T H E R M A L M O D E L S A N D C O R R E S P O N D I N G G A S P RO D U C T I O N Model A
The simplest assumption is that exposed water ice sublimates directly from the surface. The ice must be mixed with dust, considering the low albedo of cometary nuclei. If heat conduction is neglected energy conservation leads to where I eff is the effective irradiation, A V the bolometric Bond albedo, is the emissivity, σ is the Stefan-Boltzmann constant, T the surface temperature, Z the sublimation rate, and H the latent sublimation heat. The sublimation rate is given by the Hertz-Knudsen formula. Facets of a non-convex nucleus shape can be shaded but also receive energy from nearby facets. The effective irradiation includes the direct solar irradiation and the reflected visible and IR contributions, as well as the emitted temperature-dependent IR input from the facets within the field of view of the reference facet. The IR emitted and reflected components of the self-illumination are computed iteratively. This non-linear algebraic equation is solved with respect to T, which makes it possible to calculate the corresponding local gas production.
The results of computer simulations for various variants of Model A as well as approximating power functions are shown in Fig. A1.

Models B and C
Models B and C differ from Model A in that we take into account explicitly the existence of a porous dust layer on the surface and a heat flux inward due to the thermal conductivity of the material. In these models, the equation of the energy balance on the surface is replaced by a system of two equations, each of which describes the energy balance on the surface (x = 0) and at the ice/dust interface (x = L), respectively. ( A more detailed description of the models can be found in Keller et al. (2015b) and Skorov et al. (2016). It is important to note that in these models radiation thermal conductivity is taken into account. It is proportional to the pore size between the particles and the temperature to the third power. Figure A1. Total gas production as a function of the heliocentric distance R H . Simulation results are shown for: 67P SHAP7 Model A (black curve), 67P SHAP7 Model A for a zero-obliquity with (blue curve), and without self-illumination (red curve); sphere Model A (green curve). The exponent of the approximating function and the correlation coefficient is shown in the legend. Figure A2. Total gas production as a function of the heliocentric distance R H . Simulation results are shown by solid (pre-perihelion) and dashed (postperihelion) curves. Model B for SHAP7 -black colour and model C for SHAP7 -red colour. Power-law fittings are shown by crosses. The exponent of the approximating function and the correlation coefficient is given in the legends.
The results of computer simulations for Model B and Model C as well as approximating power functions are shown in Fig. A2.

Model from Hu et al. (2017)
The model presented in Hu et al. (2017) describes non-stationary heat transfer in a two-layer porous medium. This is its main difference from Models B and C, which is built on the assumption that the surface layer is in a quasi-stationary state. A detailed description of the model can be found in the cited paper. Unfortunately, due to obvious computational limitations, the authors assume the thermal Figure A3. Total gas production (molecules per sec) as a function of the heliocentric distance R H . Solid curves are from Hu et al. (2017), dashed curves are power fitting functions. Layer thickness (top to bottom): 5, 10, and 15 mm. The particle size is 1 mm. Icy factor is 0.1. The exponent of the approximating function and the correlation coefficient is shown in the legends.
conductivity to be constant. This simplification makes us cautious about the results obtained for millimetre sized particles, when radiation transport will play an important role.
Despite this remark, the results obtained for thick layers of various thicknesses are of great interest. It was in this model that cases, where the layer thickness was comparable to the depth of the diurnally varying thermal heating, were first considered. The results of computer simulations as well as approximating power functions are shown in Fig. A3. It is seen that an increase in the layer thickness leads to a noticeable increase in the exponent of the approximating function.

Model from Blum et al. (2017)
The model presented in Blum et al. (2017) takes into account radiation thermal conductivity. In order to manifest the effects associated with this, the authors consider a medium consisting of Figure A4. Total water production (molecules per sec) as a function of heliocentric distance. Solid curves are from fig. 8 in Blum et al. (2017). Thickness of the layer 1 cm, the particle size is 1 cm. The 'areal waterice coverage' is 1 (red curve), 0.1 (blue curve), and 0.05 (black curve). The power fitting functions are shown with crosses. Fitting parameters are shown in the legends.
centimetre-sized particles. We have already paid attention to the fact that the application of a continuous approach to the description of such a porous medium seems doubtful (Skorov et al. 2016(Skorov et al. , 2017. Putting aside the discussion of this issue, we note that the results for gas production shown in the article were obtained for an extreme case: the layer thickness is equal to the particle size. Obviously, in this case, for high porosity, one cannot accurately describe the boundary of the ice free layer. However, the presented results are of interest for our stu[HUK8]dy. The results of computer simulations as well as approximating power functions are shown in Fig. A4. This paper has been typeset from a T E X/L A T E X file prepared by the author.