Mathematical Simulation of Honeycomb Weathering via Moisture Transport and Salt Deposition

: Honeycomb weathering is a common phenomenon found on various rock surfaces all around the world. However, honeycomb formation mechanisms are still poorly understood. In this study, we propose a model describing moisture transport within the sandstone and erosion resulting from salt deposition during evaporation of moisture off the rock surface. The moisture transport model is based on the non-linear diffusion equation, where the volumetric moisture content is a combined parameter accounting for the moisture and gas (vapor) content. The moisture transport model accounts for the several-orders-of-magnitude decrease in moisture diffusivity, observed during drying. It was assumed that erosion occurs when the evaporation front is located close to the rock surface. The depth of erosion is proportional to the moisture ﬂow rate through the drying surface. The ABAQUS ﬁnite-element software suite was used for numerical solution of the non-linear diffusion equation. The iterative scheme of erosion simulation for different drying cycles was implemented using the Python programming language. Computations were conducted in the 2D setting for the square model with dimensions of 50 mm × 50 mm. Simulation results demonstrate the possibility of obtaining various landform shapes (honeycombs, tafoni) by varying only the value of the distribution of moisture content at the bottom side, simulating the rate of internal wetting of rock.

In spite of intensive studies conducted over a last hundred years, honeycomb formation mechanisms still remain one of the most puzzling and confusing phenomena in geology [3]. Various processes were proposed to explain honeycomb origin, such as chemical weathering, frost weathering, thermal shock, biological weathering, case hardening, salt weathering [2,[20][21][22][23][24]. Many researchers consider salt weathering the most probable cause of honeycomb formation; however, the actual mechanisms are still poorly understood. Why are thin lips formed between pits preserved [21]? Why do similar structures form on different rock types under considerably different environmental conditions [25]? Why, on the other hand, can we observe different structures, sometimes at a different stage of honeycomb formation, present on the same rock [8]? Why do honeycombs tend to form in the areas with perfect balance between wetting and drying [26]? In spite of intensive studies conducted over a last hundred years, honeycomb formation mechanisms still remain one of the most puzzling and confusing phenomena in geology [3]. Various processes were proposed to explain honeycomb origin, such as chemical weathering, frost weathering, thermal shock, biological weathering, case hardening, salt weathering [2,[20][21][22][23][24]. Many researchers consider salt weathering the most probable cause of honeycomb formation; however, the actual mechanisms are still poorly understood. Why are thin lips formed between pits preserved [21]? Why do similar structures form on different rock types under considerably different environmental conditions [25]? Why, on the other hand, can we observe different structures, sometimes at a different stage of honeycomb formation, present on the same rock [8]? Why do honeycombs tend to form in the areas with perfect balance between wetting and drying [26]?
Finding the mechanisms of honeycomb formation is a complex problem due to a very long duration of occurring processes and scale issues [27]. For example, sandstone erosion rate constitutes several dozen millimeters per 1000 years [23]. There are successful laboratory experiments substantiating the formation of spectacular sandstone landforms such as arches, alcoves, pedestal rocks, pillars [28,29], and arcades [5] as a result of erosion. However, successful laboratory experiments illustrating the whole cycle of honeycomb formation are yet to be conducted. Computer simulation may be of a help in substantiating the erosion mechanisms of sandstone landform formation. Thus, several studies simulated the formation of arches, pedestal rocks, pillars, and arcades [30][31][32][33] based on the principal idea of negative feedback between stress and erosion. In [34], the growth of tafoni was modeled by salt weathering during a sequence of wetting/drying cycles. However, no studies have yet been conducted on simulation of honeycomb formation.
In this study, the authors propose the description of various stages of honeycomb development, based on the mathematical model of moisture transport within rock and of moisture evaporation off the rock surfaces, where erosion occurs as the result of salt deposition.

Methods
The proposed physical model describes moisture transport in the rock, and is based on the non-linear diffusion equation, where the volumetric moisture content is a combined parameter accounting for the moisture and gas (vapor) content; and ( ) is the moisture diffusivity depending on the actual moisture content [35][36][37][38]: Finding the mechanisms of honeycomb formation is a complex problem due to a very long duration of occurring processes and scale issues [27]. For example, sandstone erosion rate constitutes several dozen millimeters per 1000 years [23]. There are successful laboratory experiments substantiating the formation of spectacular sandstone landforms such as arches, alcoves, pedestal rocks, pillars [28,29], and arcades [5] as a result of erosion. However, successful laboratory experiments illustrating the whole cycle of honeycomb formation are yet to be conducted. Computer simulation may be of a help in substantiating the erosion mechanisms of sandstone landform formation. Thus, several studies simulated the formation of arches, pedestal rocks, pillars, and arcades [30][31][32][33] based on the principal idea of negative feedback between stress and erosion. In [34], the growth of tafoni was modeled by salt weathering during a sequence of wetting/drying cycles. However, no studies have yet been conducted on simulation of honeycomb formation.
In this study, the authors propose the description of various stages of honeycomb development, based on the mathematical model of moisture transport within rock and of moisture evaporation off the rock surfaces, where erosion occurs as the result of salt deposition.

Methods
The proposed physical model describes moisture transport in the rock, and is based on the non-linear diffusion equation, where the volumetric moisture content θ is a combined parameter accounting for the moisture and gas (vapor) content; and D(θ) is the moisture diffusivity depending on the actual moisture content θ [35][36][37][38]: where t is the time, and ∇ is the gradient. The moisture flow model accounts for the severalorders-of-magnitude decrease in moisture diffusivity, observed during drying. In addition, as a function of moisture content, the moisture diffusivity has a clear and deep minimum at θ m moisture content, which is commonly considered to be the marker of the evaporation front location. The intensity of moisture evaporation off the rock surface decreases with the drying and goes to zero at the moisture content of θ ∞ , equal to the humidity of the air at the external surface. To describe moisture evaporation, the flux condition at the rock surface Γ is set as follows: where β is the mass transfer coefficient, h a is the relative humidity of the blown air, h m (θ) is the relative humidity of the surface material, and n is the normal to the surface Γ. Typical values of model parameters are given in Table 1. In this study, for D(θ) and h m (θ), we used linear approximation between values. Figure 2 shows the plot for h m (θ). Table 1. Typical parameters values [37].

Symbol
Value Unit humidity of the air at the external surface. To describe moisture evaporation, the flux condition at the rock surface Γ is set as follows: where is the mass transfer coefficient, ℎ is the relative humidity of the blown air, ℎ ( ) is the relative humidity of the surface material, and is the normal to the surface Γ. Typical values of model parameters are given in Table 1. In this study, for D(θ) and ℎ ( ), we used linear approximation between values. Figure 2 shows the plot for ℎ ( ).

Symbol
Value Unit 0.01 - The ABAQUS finite-element software suite (version 6.14) was used for numerical solution of the non-linear diffusion Equation (1) [39]. ABAQUS allows a researcher to obtain realistic solutions for non-linear mechanical problems quickly [40][41][42][43][44][45]. In addition, its user subroutines mechanism makes it possible to locally change the material properties and distribution of loads during computations. In this study, the USDFLD, DFLUX, The ABAQUS finite-element software suite (version 6.14) was used for numerical solution of the non-linear diffusion Equation (1) [39]. ABAQUS allows a researcher to obtain realistic solutions for non-linear mechanical problems quickly [40][41][42][43][44][45]. In addition, its user subroutines mechanism makes it possible to locally change the material properties and distribution of loads during computations. In this study, the USDFLD, DFLUX, and UEXTERNALDB user subroutines were used: USDFLD-to account for changes in material properties in finite elements during the erosion process; DFLUX-to assign the flux condition (2) at each iteration step; and UEXTERNALDB-to assign algorithm parameters before modeling and to register results of computing the moisture flow rate through the drying surface after each iteration. The user subroutines were implemented in Fortran programming language.
All computations were conducted in the 2D setting for the square model with dimensions of 50 mm × 50 mm (see Figure 3). To imitate moisture content inside the rock, the constant value for moisture content, θ 0 , is assigned at the bottom surface. For symmetry imitation, the no-penetration condition is set at the side surfaces. To imitate moisture evaporation at the rock surface, the flux conditions are set at the upper side of the model. The DC3D8 linear hexahedral elements were used in all computations discussed here. mensions of 50 mm × 50 mm (see Figure 3). To imitate moisture content inside t the constant value for moisture content, 0 , is assigned at the bottom surface. F metry imitation, the no-penetration condition is set at the side surfaces. To moisture evaporation at the rock surface, the flux conditions are set at the uppe the model. The DC3D8 linear hexahedral elements were used in all computati cussed here. Here, the sandstone erosion caused by salt deposition is considered to be th tial mechanism influencing formation of honeycomb [4,34,[46][47][48][49][50]. The process erosion is only affected by salt depositions leading to formation of subflorescence Authors believe that erosion occurs when the evaporation front is located clos rock surface:   Here, the sandstone erosion caused by salt deposition is considered to be the essential mechanism influencing formation of honeycomb [4,34,[46][47][48][49][50]. The process of rock erosion is only affected by salt depositions leading to formation of subflorescence [51,52]. Authors believe that erosion occurs when the evaporation front is located close to the rock surface: At the same time, the depth of erosion is proportional to the moisture flow rate through the drying surface q. For computations, repeating cycles of air humidity changes are modeled by changing the parameter of the relative humidity, h a . The following parameter values were considered: h a1 = 0.78; h a2 = 0.76; h a3 = 0.74. Figure 4 shows the modeling algorithm for rock surface erosion. The upper part of the model (see Figure 3) simulating the rock surface is built with the polyline consisting of line segments. When building the model, the size of each segment should be at least 0.3 mm and should not exceed 1.0 mm. After the finite element mesh is built automatically, the boundary conditions are set, including flux conditions with the current values for parameters h a1 and h m (θ). A steady state solution is found when solving the diffusion equation. The values of the relative rate of evaporation are computed for each segment, using the following equation ). Then, for each segment, we determine its displacement into the rock region as the result of erosion. The displacement is proportional to the magnitude of q. The height of the model is maintained constant and equal to 50. If necessary, the upper surface is shifted upwards. The cycle consists of k calculations for various combinations of the parameter, h a1 . . . h ak . The cycle repeats the specified number of times, N. The scheme described here is implemented in Python programming language.
Geosciences 2023, 13, 161 5 of 16 to 50. If necessary, the upper surface is shifted upwards. The cycle consists of k calculations for various combinations of the parameter, ℎ 1 … ℎ . The cycle repeats the specified number of times, N. The scheme described here is implemented in Python programming language.

Results and Discussion
In order to analyze the influence of boundary conditions and rock surface morphology on moisture transport, position of the evaporation front, and the rate of evaporation (see Figure 5), simulations were performed with the rectangular 2D object. At the bottom side, we set conditions of constant distribution of moisture content 0 to simulate the rate of internal wetting of rock. At the upper side we set flux condition to describe moisture evaporation off the rock surface. At the left and right sides, the conditions of zero moisture transport are set. During simulations, we analyzed the data on hydraulic field distribution obtained after reaching the steady state condition. First, we consider the rectilinear upper edge simulating the smooth surface of the rock (see Figure 5a, Supplementary Figure S1) at different levels of rock saturation, i.e., the moisture content 0 (Supplementary Figure S1). The results show that the rate of moisture evaporation off the rock surface, q, is maximal for 0 > 0.118 and, in certain cases, constitutes 1.0. Within the range of 0.118 > 0 > 0.1, the rate of evaporation, q, decreases sharply and reaches 0.37, and then, at 0.1 > 0 > ∞ , a smooth fall-off to zero can be observed (see Figure 5b). The evaporation front appears at 0 < 0.118; the distance between the evaporation front and the upper edge increases with a decrease in 0 ; i.e., the depth of the vapor zone increases (Supplementary Figure

Results and Discussion
In order to analyze the influence of boundary conditions and rock surface morphology on moisture transport, position of the evaporation front, and the rate of evaporation (see Figure 5), simulations were performed with the rectangular 2D object. At the bottom side, we set conditions of constant distribution of moisture content θ 0 to simulate the rate of internal wetting of rock. At the upper side we set flux condition to describe moisture evaporation off the rock surface. At the left and right sides, the conditions of zero moisture transport are set. During simulations, we analyzed the data on hydraulic field distribution obtained after reaching the steady state condition. First, we consider the rectilinear upper edge simulating the smooth surface of the rock (see Figure 5a, Supplementary Figure S1) at different levels of rock saturation, i.e., the moisture content θ 0 (Supplementary Figure S1). The results show that the rate of moisture evaporation off the rock surface, q, is maximal for θ 0 > 0.118 and, in certain cases, constitutes 1.0. Within the range of 0.118 > θ 0 > 0.1, the rate of evaporation, q, decreases sharply and reaches 0.37, and then, at 0.1 > θ 0 > θ ∞ , a smooth fall-off to zero can be observed (see Figure 5b). The evaporation front appears at θ 0 < 0.118; the distance between the evaporation front and the upper edge increases with a decrease in θ 0 ; i.e., the depth of the vapor zone increases (Supplementary Figure S1  In order to model honeycomb weathering during erosion in accordance with the proposed method, the square computational model is used, where the upper side imitates the surface of rock and has a random structure with 10 mm-wide and 4 mm-deep pits (see Figure 3). During simulations, only the value of constant distribution of moisture content 0 at the bottom side was changed, simulating the rate of internal wetting of rock. Figure 6 shows simulation results for constant distribution of moisture content at the bottom side set to 0 = 0.115. Moisture content over the whole volume exceeds , showing that the model is fully wet, and no erosion takes place (3). Figure 7 shows simulation results after 40 iterations with the constant distribution of moisture content at the bottom side set to 0 = 0.110. A virtually flat surface can be observed here. This can be explained that only protruding parts dry out below , and hence, only protruding parts undergo erosion. The evolution of a pitted surface into a virtually flat one is shown in Supplementary Video S1. In order to model honeycomb weathering during erosion in accordance with the proposed method, the square computational model is used, where the upper side imitates the surface of rock and has a random structure with 10 mm-wide and 4 mm-deep pits (see Figure 3). During simulations, only the value of constant distribution of moisture content θ 0 at the bottom side was changed, simulating the rate of internal wetting of rock. Figure 6 shows simulation results for constant distribution of moisture content at the bottom side set to θ 0 = 0.115. Moisture content over the whole volume exceeds θ m , showing that the model is fully wet, and no erosion takes place (3). Figure 7 shows simulation results after 40 iterations with the constant distribution of moisture content at the bottom side set to θ 0 = 0.110. A virtually flat surface can be observed here. This can be explained that only protruding parts dry out below θ m , and hence, only protruding parts undergo erosion. The evolution of a pitted surface into a virtually flat one is shown in Supplementary Video S1.   Figure 8 shows simulation results after 40 iterations; the constant distribution of moisture content at the bottom side is set to 0 = 0.109. The formation of the isolated single lip can be observed where the protruding part undergoes significant drying, resulting in a considerable reduction in the evaporation rate. Simultaneously, a more intense evaporation takes place at other regions. Thus, the lip forms from the initially protruding part that undergoes significant drying as compared to other rock surface. The formation of the lip corresponds to the formation of honeycombs. The evolution of the initial surface with pits into a virtually flat surface with a single isolated lip is shown in Supplementary Video S2. We can see that the lip stays dry all the time, as the evaporation front is located at the foundation of the lip. Thus, no erosion occurs, and the lip retains its shape. At the lesser value of 0 = 0.108 the number of lips increases due to a larger number of protruding parts that undergo drying, resulting in a significant reduction in   Figure 8 shows simulation results after 40 iterations; the constant distribution of moisture content at the bottom side is set to 0 = 0.109. The formation of the isolated single lip can be observed where the protruding part undergoes significant drying, resulting in a considerable reduction in the evaporation rate. Simultaneously, a more intense evaporation takes place at other regions. Thus, the lip forms from the initially protruding part that undergoes significant drying as compared to other rock surface. The formation of the lip corresponds to the formation of honeycombs. The evolution of the initial surface with pits into a virtually flat surface with a single isolated lip is shown in Supplementary Video S2. We can see that the lip stays dry all the time, as the evaporation front is located at the foundation of the lip. Thus, no erosion occurs, and the lip retains its shape. At the lesser value of 0 = 0.108 the number of lips increases due to a larger number of protruding parts that undergo drying, resulting in a significant reduction in  Figure 8 shows simulation results after 40 iterations; the constant distribution of moisture content at the bottom side is set to θ 0 = 0.109. The formation of the isolated single lip can be observed where the protruding part undergoes significant drying, resulting in a considerable reduction in the evaporation rate. Simultaneously, a more intense evaporation takes place at other regions. Thus, the lip forms from the initially protruding part that undergoes significant drying as compared to other rock surface. The formation of the lip corresponds to the formation of honeycombs. The evolution of the initial surface with pits into a virtually flat surface with a single isolated lip is shown in Supplementary Video S2. We can see that the lip stays dry all the time, as the evaporation front is located at the foundation of the lip. Thus, no erosion occurs, and the lip retains its shape. At the lesser value of θ 0 = 0.108 the number of lips increases due to a larger number of protruding parts that undergo drying, resulting in a significant reduction in the evaporation rate (see Figure 9). The evolution of the initial surface with pits into a virtually flat surface with several isolated lips is shown in Supplementary Video S3. the evaporation rate (see Figure 9). The evolution of the initial surface with pits into a virtually flat surface with several isolated lips is shown in Supplementary Video S3.   Figure 10 shows simulation results after 40 iterations; constant distribution of moisture content at the bottom side is set to 0 = 0.106. We can observe the increasing number and thickness of lips. The evolution of the initial surface into a virtually flat one with several isolated thicker lips is shown in Supplementary Video S4. Further reduction in the value of constant distribution of moisture content at the bottom side to 0 = 0.100 results in a growing thickness of lips and merging of some adjacent lips (see Figure 11 and Supplementary Video S5). the evaporation rate (see Figure 9). The evolution of the initial surface with pits into a virtually flat surface with several isolated lips is shown in Supplementary Video S3.   Figure 10 shows simulation results after 40 iterations; constant distribution of moisture content at the bottom side is set to 0 = 0.106. We can observe the increasing number and thickness of lips. The evolution of the initial surface into a virtually flat one with several isolated thicker lips is shown in Supplementary Video S4. Further reduction in the value of constant distribution of moisture content at the bottom side to 0 = 0.100 results in a growing thickness of lips and merging of some adjacent lips (see Figure 11 and Supplementary Video S5).     Figure 12 shows simulation results after 67 iterations; the constant distribution of moisture content at the bottom side is set to 0 = 0.06. The evaporation front is located closer to the bottom of the pit than in other regions. That is why the erosion actively develops at the bottom of the pit. Here, one can observe the formation of the deep isolated pit (tafoni) from the initial pit. In other regions, the erosion process is less pronounced because the evaporation front is removed from the surface. The formation of the deep isolated pit (tafoni) from the initial pit is shown in Supplementary Video S6. At 0 = 0.03, the evaporation front is located at a significant distance from the surface (see Figure 13), and no erosion takes place (3).   Figure 12 shows simulation results after 67 iterations; the constant distribution of moisture content at the bottom side is set to 0 = 0.06. The evaporation front is located closer to the bottom of the pit than in other regions. That is why the erosion actively develops at the bottom of the pit. Here, one can observe the formation of the deep isolated pit (tafoni) from the initial pit. In other regions, the erosion process is less pronounced because the evaporation front is removed from the surface. The formation of the deep isolated pit (tafoni) from the initial pit is shown in Supplementary Video S6. At 0 = 0.03, the evaporation front is located at a significant distance from the surface (see Figure 13), and no erosion takes place (3).  Figure 12 shows simulation results after 67 iterations; the constant distribution of moisture content at the bottom side is set to θ 0 = 0.06. The evaporation front is located closer to the bottom of the pit than in other regions. That is why the erosion actively develops at the bottom of the pit. Here, one can observe the formation of the deep isolated pit (tafoni) from the initial pit. In other regions, the erosion process is less pronounced because the evaporation front is removed from the surface. The formation of the deep isolated pit (tafoni) from the initial pit is shown in Supplementary Video S6. At θ 0 = 0.03, the evaporation front is located at a significant distance from the surface (see Figure 13), and no erosion takes place (3).     The simulation results shown in Figures 6-13 and in Table 2 Table 2   To illustrate the simulation results, Figure 14 shows the photographs of the edge of the rock located in Apolena Rock City, CZ, where various landform shapes are represented. The lower right part of the rock is close to the earth and is, presumably, very wet. Hence, no erosion occurs there (see Figure 6). Higher up, the internal wetting of rock decreases. First, we can see the flat surface of rock (see Figure 7). Further on, as the height increases, we can observe the region with honeycombs (see Figure 9). Closer to the edge of the rock, the intense drying of the rock can be observed, resulting in the formation of tafoni (see Figure 12). At the very edge of the rock, the material is presumably very dry, and no erosion occurs (see Figure 13).  To illustrate the simulation results, Figure 14 shows the photographs of the edge of the rock located in Apolena Rock City, CZ, where various landform shapes are represented. The lower right part of the rock is close to the earth and is, presumably, very wet. Hence, no erosion occurs there (see Figure 6). Higher up, the internal wetting of rock decreases. First, we can see the flat surface of rock (see Figure 7). Further on, as the height increases, we can observe the region with honeycombs (see Figure 9). Closer to the edge of the rock, the intense drying of the rock can be observed, resulting in the formation of tafoni (see Figure 12). At the very edge of the rock, the material is presumably very dry, and no erosion occurs (see Figure 13). The results presented in this study are the first attempt to describe the process of honeycomb formation via numerical simulation. The results obtained here can provide answers to some questions formulated in the Introduction section. Thin lips formed between pits are retained because the evaporation front is located at the foundation of lips, and no intensive evaporation takes place inside the pits. This phenomenon is not inherent in specific types of rocks, but is determined by the physical processes of liquids evaporation. Taking into account that the simulations were conducted with typical properties of the rock and its environment, we can expect that the formation of structures on different The results presented in this study are the first attempt to describe the process of honeycomb formation via numerical simulation. The results obtained here can provide answers to some questions formulated in the Introduction section. Thin lips formed between pits are retained because the evaporation front is located at the foundation of lips, and no intensive evaporation takes place inside the pits. This phenomenon is not inherent in specific types of rocks, but is determined by the physical processes of liquids evaporation. Taking into account that the simulations were conducted with typical properties of the rock and its environment, we can expect that the formation of structures on different rock types takes place under considerably different environmental conditions. The formation of honeycombs is related to the following physical phenomena: decrease in moisture diffusivity of the rock by several orders of magnitude under drying and decrease in evaporation rate during drying. Rock erosion is only affected by salt depositions leading to the formation of subflorescence that occurs when the evaporation front is close to the surface. That is why the formation of honeycombs can be expected in the presence of physical phenomena described earlier and under specific moisture content in the rock. We also proposed the example describing the formation of various shapes at the same rock, where the differences in surface structures are explained by the non-uniform distribution of moisture inside the rock. Simulation results demonstrate that honeycombs tend to form in the areas having perfect balance between wetting and drying.
Additionally, it is necessary to note and comment on certain features of the proposed numerical algorithm. There is no specific time in the calculations. We can say that equal time intervals are allocated for each case. The distinctive feature of this simulation is that we assume a constant value of the moisture content on the lower part of the model (the inner part of the rock). During computations of each case, a steady state solution is sought. In real time, the steady state can be achieved in 5-20 h, depending on the air humidity. Each case can be interpreted as a season lasting several months; i.e., different time scales are used in computations. As there is a separate calculation for each case, the humidity changes occur instantly. The computations are based on the assumption of a stationary condition for the moisture content on the lower surface. Therefore, after reaching a stationary solution, the position of the drying front does not change. The amount of subflorescence precipitated is calculated simply by multiplying local distribution of the moisture flow rate over a period of time. Thus, the variable length of wetting and drying periods can be easily taken into account.
The results presented here are qualitative. In order to quantitatively describe the formation of various shapes for concurrent rocks and to predict erosion in, for example, structures of building materials [53,54], it is necessary to measure material properties, such as moisture diffusivity D(θ) and the relative humidity of the surface material h m (θ), as functions of moisture content θ [55][56][57][58]. It is also necessary to validate the proposed physical and numerical model based on specific objects [59], by measuring seasonal changes in evaporation rate [60,61], the distribution of moisture and salts in the rock [62], and the position of the evaporation front [63]. It is also necessary to identify the influence of other factors that affect rock decay (chemical weathering [23], frost weathering [64], thermal shock [65], biological weathering [8], case hardening [4], erosion caused by the wind pressure [21], etc.). Additionally, the decay rate model needs to be improved by taking into account the petrography, mineralogy, and structural heterogeneity of the stone. The authors also intend to adapt the proposed model for the simulation of 3D structures, since this work does not describe the modeling of the formation of a regular honeycomb structure. These questions will be addressed in further studies.
This mathematical model can be applied in the forecasting of salt weathering of porous building materials, prediction of the life span of building structures, and developing methods for the conservation of historical monuments and our cultural heritage [53,66,67]. In addition, the developed algorithm for solving the non-linear diffusion equation with the evolution of the solution domain and boundary conditions can be applied to a wide range of natural and biological processes including neural networks, vascular networks, slime mold, plant routes, fungi mycelium, bone remodeling, bacterial biofilms [68] and colony patterns [69,70], the formation of the skeletons of microorganisms [71,72], fulgurites [73], and snowflakes [74], etc.

Conclusions
The proposed model demonstrates the formation of various landform shapes (honeycombs, tafoni) through the variation of moisture content, simulating the rate of internal wetting of the rock. At high values of moisture content in the rock, the evaporation front is located close to the surface, and thus, almost the entire surface remains wet, and erosion occurs only on the protruding parts of the surface. This results in smoothing of the surface. However, at low moisture content, the evaporation front moves away from the surface, and a more intense flow is observed at the bottom of the pits, leading to their deepening. Thus, at different values of moisture content inside the rock, two opposite phenomena can be observed-smoothing of the surface and deepening of the pits. At intermediate values of moisture content, the interaction of these opposite phenomena results, accordingly, in the formation of lips and honeycombs.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.