Parametric study of lotus-type pore shape in solid subject to Henry's laws at interfaces

Mechanisms of the length and maximum radius of lotus-type or single pores in ice or nonmetals satisfied by Henry's law at gas-liquid interfaces dissolved by a gas during unidirectional solidification are rigorously investigated and supported by a Table from algebraic predictions involving different dimensionless working parameters. Lotus-type porous materials characterized by directional properties have been often used as functional materials in food, biomedical, and micro- and nano-technologies. Following previous work taking into account solute amount and transport within the pore, and concentration boundary layers on the advancing solid-liquid interface and bubble cap, and the Young-Laplace equation and Henry's law at liquid-gas interfaces, the algebraic study further provides a Table for a quantitative and extensive understanding of different mechanisms of length and maximum radius. Dimensionless parameters include solute transport parameters of Henry's law constant, mass transfer coefficient, partition coefficient, solute gas amount in imposed ambient, and solute transport parameter, and fluid and thermal parameters of solidification rate, imposed gas pressure, hydrostatic pressure, and geometrical parameter of inter-pore spacing. The controlling of the shapes of lotus-type pores is achieved by a good comparison between predicted maximum diameter and inter-pore spacing during freezing of water dissolved by oxygen gas.


Introduction
Materials containing gas porosity such as ice or nonmetals are encountered in various fields. For example, ice containing gases such as nitrogen and oxygen has been often used in food and fishery and preservation industries [1]. In view of strong oxidation, ice containing ozone has widely been used for deodorization, sterilization, purification of water, and prolongation of shelf-life foods [2,3]. Ozone gas can be dissolved into water with an ozone injection system, and fed into an ice making machine that will freeze water to produce ozonated ice. More than eighty years ago, research in the commercial fish industry was conducted in France that showed a 33% extension in shelf-life of fresh fish stored on ozonated ice when compared with ice produced from regular water in the holds of fishing vessels. Since gas amount in the ice is limited by gas solubility in water, Inada et al. [1] thus proposed a microbubble method to enhance ice containing gas pores. Microbubbles of a desired gas are dispersed, and then obtained images of the structural features of pores formed in the ice. The air amount in ice prepared by the microbubble method was three times higher than the solubility of air in water prepared by the conventional method. Matsumoto et al. [4] also studied pore formation during freezing water containing ozone micro-bubbles finding the effects of surfactant, inclination angles of the cooling plate on characteristics of ice containing oxygen micro-bubbles.
Ice containing pores filled with methane, oxygen or carbon dioxide has been often seen and well treated in geophysical sciences [5][6][7], and fundamental sciences [8,9]. In view of reflecting sunlight and obstructing thermal energy, growth and decay of ice in ponds, lakes and sea ice, which contain air pores are of interest. The occurrence of bubbles in fresh water ice relates to the gas amount in the pond/lake water prior to freezing over the surface. As freezing progresses, the interface concentration of dissolved gases surpasses a critical value, the water at the interface becomes supersaturated, and gas bubbles nucleate and grow to a visible size. When the sea water starts to freeze, most of salt amount is discharged from ice crystals into liquid below because the solubility limit of salt in the crystal lattice of ice is very small. However, a number of small volumes of saline water, called brine regions, and air pores remain in crystal grains of ice. Sea ice is thus also a porous material which is composed of pure ice crystals, brines and air inclusions. Crabeck et al. [6] discussed that formation of sea ice in winter can trap salts and CO 2 in brine from the ocean. CO 2 trapped with salt is transformed into solid rock. CO 2 also trapped in bubbles can rise from the bottom of the sea ice to the surface through the brine or air channels. Through the rising bubbles and sinking brine, the sea ice loses a lot of CO 2 . As a result, when the sun comes back in the spring, the sea ice no longer holds as much CO 2 .Sea ice can again absorb lots of CO 2 from the ambient. Sea ice thus can help the ocean to absorb CO 2 to fight climate change.
Shape development of the lotus-type pores in ice during upward freezing of water-carbon dioxide solutions were observed by Murakami and Nakajima [8]. Columnar pores were formed at low growth rates. The length of the columnar pores became shorter as the growth rate increased or the carbon dioxide concentration increased. The length of the columnar pores and periodic formation of pores were explained in terms of the rate of supply of carbon dioxide to the growing pores from the surrounding liquid. Yoshimura et al. [9] observed that lotus-type pores during solidification of water containing oxygen gas changed from egg, cylindrical into bifurcated cylindrical shapes as solidification rate increased. Theoretical analysis was conducted to scale solute transport in radial direction within the concentration boundary layer ahead of the advancing solid-liquid interface into pore, satisfied by Henry's law applicable to liquid containing a molecular solute gas with small concentration at the bubble cap. Diameter and inter-pore spacing were found to be inversely proportional product of solidification rate with square root of ambient pressure, agreeing with measurements. Other than proposing a universal phase diagram and simultaneous first-order ordinary differential equations [10], Wei and Huang [11] proposed a polynomial with three degrees to predict isolated pore shapes in solid by accounting for either solute transport from the pore across the bubble cap into surrounding liquid or that from the surrounding liquid into the pore in the early stage, denoted by Cases 1 and 2, respectively. Recently, Wei and Lee [12] combined Cases 1 and 2 interrelated by interaction angle and solute transport parameter to algebraically predict morphology of lotus-type pores. The predicted pore shapes were in good agreements with available experimental results during freezing of water filled with oxygen gas and carbon dioxide, and liquid copper imposed by hydrogen gas.
In this study, morphology of lotus-type pores satisfied by Young-Laplace equation and Henry's law at the bubble cap and top free surface during unidirectional solidification are parametrically presented. Instead of solving simultaneous first-order differential equations [10,11], this study extends previous work [12] to provide a table to get deep insight into understanding of dimensionless working parameters responsible for different shapes of lotus-type pores.

System model and analysis
Lotus-type pores with transverse and longitudinal cross-sections in maximum radius, r B90 , and inter-pore spacing, w, are illustrated in Fig. 1(a) and (b), respectively. Since solute rejected by the advancing solid-liquid interface results in supersaturation for bubble nucleation, concentration boundary layers occur along the advancing solid-liquid interface and on the bubble cap. The major assumptions made are as follows: 1. The pore grows unidirectionally. Furthermore, the variations of transport variables and pore shape in the circumferential direction are much smaller than those in the axial direction. This can be confirmed by the lotus-type pores in cylindrical shape (see Fig. 1(a)).
2. The tiny bubble with spherical caps at the top and bottom is characterized by a small Bond number. Bond number Bo ≡ρ l gR 2 0 / σ representing the ratio between hydrostatic pressure and capillary pressure can be as small as 10 − 5 in most cases, based on a typical bubble radius R 0 = 10 − 5 m. 3. The tiny pore is in a lumped system. 4. Fluid flow is dominated by force convection rather than free convection driven by heat transfer or mass transfer [13]. Velocity of free convection driven by heat transfer in the thermal boundary layer can be scaled by considering balance between viscous force and free convection due to heat transfer, μv T,NC /δ 2 T ∼ ρβ T gΔT, and thermal convection and conduction, v T,NC /H ∼ α l /δ 2 T . Velocity of free convection driven by heat transfer thus yields [14] v T,NC ∼ where Rayleigh number due to temperature Ra T ≡ ρβ T gΔTH 3 /α l μ ∼ 10 − 5 based on temperature expansion coefficient β T ~5× 10 − 5 1/ K, dynamic viscosity μ ∼ 1.7 × 10 − 3 kg/m − s, and thermal diffusivity α l ∼ 1.3 × 10 − 7 m 2 /s of water [14], length of the physical domain H ∼ 10 − 5 m and temperature difference ΔT ~1 K. Similarly, considering balance between viscous stress and free convection due to concentration μv C,NC /δ 2 C ∼ ρβ C gΔC and balance between convection and diffusion v C,NC /H ∼ D l /δ 2 C in the concentration boundary layer gives velocity of free convection driven by mass transfer [14] v C,NC ∼ D l H Ra where Rayleigh number due to concentration Ra C ≡ ρβ C gΔCH 3 /μD l ∼ 10 − 5 , based on concentration expansion coefficient of water due to oxygen gas or carbon dioxide β C ∼ 10 − 4 m 3 /kg, solute diffusivity D l ∼ 10 − 9 m 2 /s [14], and difference in solute concentration ΔC ~10 mol/m 3 . Referring to Eqs. (1) and (2) indicates that velocity of free convection driven by either heat transfer or mass transfer is much smaller than velocity of force convection greater than 0.01 m/s [13].
where Marangoni number Ma T ≡ − (dσ /dT)ΔTH/μα l ~7, Prandtl number Pr ≡ μ/ρα l ~13 [14], and surface tension coefficient due to temperature dσ/dT ~ 1.5 × 10 − 4 N/m-K [16]. Similarly, viscous stress balanced by Marangoni force due to concentration, μv C,FC / δ FC ∼ (dσ /dC)ΔC/H, and inertia force balanced by viscous stress ρv 2 where Marangoni number Ma C ≡ − (dσ /dC)ΔCH/μD ~60, surface tension coefficient due to concentration dσ/dC ∼ 10 − 6 Nm 2 / mol [17,18], and Schmidt number Sc ≡ μ/ρD ~450 [14]. Eqs. (3) and (4) indicates that Marangoni force can be conditionally considered, provided that velocity due to force convection is not much higher than 0.1 m/s. 6. The weakness of this model can be attributed to the strong effects of thermocapillary convection due to temperature on fluid flow, as stated in item 5. As a result, liquid pressure in the Young-Laplace equation deviates from hydrostatic pressure, which is only a function of depth in the liquid layer and gravitational acceleration. However, integration method used in this model to evaluate solute content in the system is acceptable, even though convection is accounted [19]. 7. Tangential and normal stresses exerted by fluid flow are thus negligible in this static system. Liquid pressure is identical to hydrostatic pressure. 8. Inter-pore spacing which is a result of constitutional supercooling is specified [20].
With the above assumption, the entrapped bubble in solid is divided into three segments, defined by initial contact angle, and inclination angle or contact angle of 90 • , as sketched by Fig. 2 [12]. The pore shape is closely related and described by the second segment.
Other degree of polynomials can be found in previous work [21].

Solute conservation
Determination of length and maximum radius of the second segment in Eqs. (5)-(11) accounts for solute conservation in a controlled system.

Total solute conservation
Solute molecules in the controlled system subject to adiabatic boundary conditions are conserved and given by where the left-hand side represents conserved total solute molecules including solute molecules within the entrapping bubble, and boundary layers on the cap and advancing solid-liquid interface, respectively. Intersection angle ϕ cr between two boundary layers is illustrated in Fig. 3. Integrals in Eq. (12) can be evaluated by substituting concentration profiles within boundary layers on the cap and advancing solid-liquid interface [12] C where concentration boundary layer thicknesses [11,12].
Pore volume in Eq. (12) is where integral can be approximated by
Henry's law at the bubble cap and top surface are, respectively [22].
Dimensionless length of the second segment by substituting Eqs. (13)-(17) into Eq. (12) gives Fig. 4. Prediction and measurement of (a) maximum diameter and (b) inter-pore spacing versus solidification rate during freezing of water dissolved by oxygen gas with different pressures [12].
where dimensionless functions and Eqs. (24) and (25) are, respectively, the Young-Laplace equation in the pore at initial time and Henry's law at the top surface.

Solute concentration in pore
Changing rate of solute amount in the entrapping bubble at contact angle of 90 • yields where terms on the right-hand stand for solute transport from boundary layers on the cap and advancing solid-liquid interface, respectively. The dimensionless maximum radius from Eq. (26) thus yields with Eq. (27) with algebraic coefficients from Eqs. (28)-(32) can be exactly and analytically obtained by fundamental mathematics.

Results and discussion
Lotus-type pores in solid are algebraically predicted for different independent dimensionless working parameters based on selected typical values of K Hc = K H∞ = 34, U 0 = 0.11, U 90 = 0.035, k p0 = 0.0088, k p90 = 0.008, h D0 = 0.09, h D90 = 0.03, p a = 790, Bo(h B0 − z B0 ) = 7.86, γ = 1, w = 2.5, ϕ B0 = 140 • , and η = 1. Instead of solving simultaneous first-order ordinary equations [10], algebraic equations explicitly gain deep insight into mechanisms of length and maximum radius of lotus-type or single pores. Henry's law constants are allowed to be different at the bubble cap and top free surface, since the latter is higher than the former due to a decrease in surface curvature [23]. Solidification rates, mass transfer coefficients and partition coefficients can be varied during solidification. Since concentration boundary layer thickness increases as a result of decrease in solidification rate, mass transfer coefficient on the bubble cap decreases in the course of solidification, as can be seen from Eq. (14). Because effective partition coefficient increases with solidification rate [24], partition coefficient at inclination angle of 90 • is specified to be less than that at initial contact angle. Partition coefficient can also be given by   Substituting Eq. (14) into Eq. (33) shows that partition coefficient is relatively constant and greater than equilibrium partition coefficient. Since partition coefficient is defined by the ratio between solute concentrations in solid and liquid at the advancing solidliquid interface, its value can be affected by convection, boundary layers, interface morphology, etc. [24].
The length of lotus-type pores can be approximated by the length of the second segment determined by Eq. (18), governed by solute in the boundary layer on the entrapping bubble at initial time and contact angle of 90 • , namely, , and solute in boundary layer on the advancing solid-liquid interface at initial time and contact angle of 90 , and F 5 (k p90 , r B90 , w, U 90 ) are, respectively, related to thicknesses of concentration boundary layers on the cap and advancing liquid-solid interface (see Eqs. (19), (20), (22) and (23)). Therefore, pore length increases as solute amount in boundary layers on the cap and advancing solid-liquid interface at initial time increases, whereas that at inclination angle of 90 • decrease. Solute in top and bottom spherical caps at contact angle of 90 • is referred to p g90 F 3 (ϕ B0 , r B90 ) from Eq. (21). The maximum radius is, however, calculated from Eq. (27) related to time change in solute amount in the entrapping bubble as a result of solute transfer rates through the cap from boundary layers at contact angle of 90 • .
Predicted maximum diameter and inter-pore spacing versus solidification rate during unidimensional freezing of water subject to different imposed oxygen gas pressures at the top free surface agree with available experimental results [9], as presented in Fig. 4(a)   Fig. 7. Predicted lotus-type pore shapes affected by dimensionless solidification rates at time corresponding to initial time. and (b), respectively. Increases in solidification speed and imposed solute gas pressure at the top surface decrease the maximum diameter and inter-pore spacing. Polynomials with three, four and five degrees to predict morphology of lotus-type pores are shown in Fig. 5. It can be seen that a polynomial with five degrees used in this work is sufficiently accurate to delineate morphology of lotus-type pores.
Dimensionless length of lotus-type pores increases with dimensionless Henry's law constant at the cap, as can be seen from Fig. 6  (a). Referring to Table 1, this is attributed to decrease in solute amount in the boundary layer on the cap at contact angle of 90 • . Difference in concentrations across the boundary layer between the cap and top free surface decreases as Henry's law constant at the cap increase. On the other hand, pore length increases as Henry's law constant at the top surface decreases, as shown in Fig. 6(b). This is because of increase in solute amount in the boundary layer on the advancing solid-liquid interface at initial time. Predicted result also agrees with that of a single isolated pore and specification of solute concentration in liquid far from the pore [10][11][12].
The effects of dimensionless solidification rate at the initial time on the size of lotus-type pores are shown in Fig. 7. An increase in dimensionless solidification rate at the initial time from 0.1 to 0.11 for solidification rate at contact angle of 90 • of 0.035 is seen to decrease pore length. This is because of decrease in solute amount in the boundary layer on the advanced solid-liquid interface at initial time (see Table 1).
The effects of dimensionless mass transfer coefficient at initial contact angle on the pore shape are shown in Fig. 8. Dimensionless  The effects of partition coefficient at initial contact angle on the shape of lotus-type pores are shown in Fig. 9. A decrease in partition coefficient at initial contact angle significantly increases the length of the pore, as a result of an increase in solute in the boundary layer on the advancing solid-liquid interface at the initial time.
The effects of mole fraction of solute gas at the top free surface on the size of lotus-type pores are shown in Fig. 10. Length of lotustype pores increases with mole fraction of solute gas on the top free surface. This is attributed to increase in solute amount in the boundary layer on the advancing solid-liquid interface at initial time (see Table 1). It is worthy of mentioning that increase in mole fraction of solute gas on the top free surface decreases solute amount in concentration boundary layer on the cap whereas enhances that on the advancing solid-liquid interface (see Table 1). Mole fraction of solute gas on the top free surface exhibits stronger effects on solute amount in the concentration boundary layer on the solid-liquid interface than that on the bubble cap. The effects of mole fraction of solute gas at the top free surface on solute amount in the concentration boundary layer on the advancing solid-liquid interface at initial time are also stronger than that at contact angle of 90 • .
An increase in dimensionless imposed pressure at the top surface increases pore length, as presented in Fig. 11. This is because of an increase in solute amount in the boundary layer on the advancing solid-liquid interface at initial time. An increase in dimensionless  hydrostatic pressure decreases dimensionless length of lotus-type pores, as shown in Fig. 12. This is attributed to an increase in solute amount in the boundary layer on the cap at contact angle of 90 • .
Figs. 13 and 14 show that pore length increases with initial contact angle and dimensionless inter-pore spacing. This is a consequence of an increase in solute amount within the boundary layer on the advanced solid-liquid interface at initial time. An increase in solute transport parameter increases radius and length of the pore, as shown in Fig. 15. In view of enhanced solute transport from the boundary layer on the advancing solid-liquid interface into the pore, increase in solute transport parameter decreases solute amount in the boundary layer on the advancing solid-liquid interface at contact angle of 90 • .

Conclusions
The conclusions drawn are as follows: 1. The lotus-type pore shapes characterized by the length and maximum radius affected by working parameters are systematically analyzed from algebraic expressions and Table provided. Prediction and available measurement of diameter and inter-pore spacing during freezing of water dissolved by oxygen gas are in good agreement.  2. Length of lotus-type pores are usually controlled by solute content in the pore and concentration boundary layer on the solidification front. Solute content in the concentration boundary layer on the cap at contact angle of 90 • is more important than that at initial time. Solute content in the concentration boundary layer on the initial time is more significant than that at contact angle of 90 • . 3. Length and maximum radius exhibit different mechanisms. The former is determined by maintaining total solute amount in the system, composed of the pore and boundary layers on the cap and advancing solid-liquid interface between initial time and time corresponding to contact angle of 90 • . The latter is calculated by considering changing rate of solute amount in the pore due to solute transport from boundary layers on the cap and advancing solid-liquid interface at contact angle of 90 • . 4. Dimensionless length of lotus-type pores increases with Henry's law constant at the cap. Referring to the Table provided, this is primarily due to decrease in solute amount within the boundary layer on the cap at contact angle of 90 • . On the other hand, length of lotus-type pores increases as Henry's law constant at the top surface decreases, resulting from increase in solute amount in the boundary layer on the advancing solid-liquid interface at initial time. 5. Dimensionless length of lotus-type pores increases as solidification rate decreases at initial time. This is because increase in thickness results in increase of solute amount in the boundary layer on the advancing solid-liquid interface at initial time. Similar results can be found for partition coefficients at initial time. 6. Dimensionless length of lotus-type pores increases as dimensionless mass transfer coefficients at initial time decreases. This is because increase in thickness increases solute amount in the concentration boundary layer on the cap at initial time. 7. Dimensionless length of lotus-type pores increases with initial contact angle, dimensionless imposed pressure and mole fraction of solute gas on the top free surface, and inter-pore spacing. This is primarily attributed to enhanced solute amount in boundary layer on the advancing liquid-solid interface at initial time. 8. Dimensionless length of lotus-type pores increases as dimensionless hydrostatic pressure decreases, as a result of decrease in solute amount in the boundary layer on the cap at contact angle of 90 • . 9. Dimensionless pore length increases with solute transport parameter, because of decrease in solute within the boundary layer on the advancing solid-liquid interface at contact angle of 90 • .
Author contribution statement P. S. Wei: Conceived and designed the experiments; Analyzed and interpreted the data; Wrote the paper. B. Y. Lee: Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data.

Data availability statement
No data was used for the research described in the article.