In restoration work, it is clear that the shallow fractures on the surface are easily repaired, while for deeper fractures, it is necessary to simulate and verify whether the restoration fluid will achieve the desired results. In the early stages of restoration work, the flow of the restoration fluid through the cracks needs to be considered. The classical flow dynamics behaviour of fluids is summarised in Navier-Stokes Equations, expressed as follows.

$\frac{\partial (\rho u)}{\partial t}+\vec{u}\cdot\nabla(\rho \vec{u})=-\nabla p+\nabla \cdot \tau +f$

u denotes fluid velocity, p denotes fluid pressure, τ denotes shear stress, f denotes the external force term and ρ denotes fluid density. For an incompressible (fluid density remains constant) Newtonian fluid (the effect of viscosity on fluid velocity is linear), the NS equation reduces to:

$\frac{\partial u}{\partial t}+\vec{u}\cdot \nabla\vec{u}=-\frac{1}{\rho}\nabla p+\mu \cdot \Delta u+\frac{f}{\rho}$

μ denotes the fluid viscosity, the magnitude of which is influenced by two factors, the concentration of ions in the solution and the volume fraction of the suspended solids. When the concentration of the solid phase in the suspension is considered alone, the relative viscosity can be expressed as follows [1]:

$\mu_{r1}=\frac{\mu}{\mu_{water}}=(\frac{1-\Phi}{1-\frac{\Phi}{\phi_{max}}})^{\frac{C_1\cdot\phi_{max}}{1-\phi_{max}}}$

When considering the concentration of ions in solution alone, the relative viscosity can be expressed as follows [2]:

$\mu_{r2}=1+A\cdot c^{\frac{1}{2}}+B\cdot c$

A and B are correlation coefficients with different values for different substances, C₁ is taken as 2.5 and Φ denotes the volume fraction of the solid phase, which for Bacillus subtilis can be calculated as follows:

$\Phi=c_N\cdot V_{Bacillus}\times 100%$

Where c_N is the concentration of Bacillus subtilis in the restoration solution, taken as 10⁸ per ml, and V_Bacillus is the volume of individual Bacillus subtilis, taken as 0.75 μm x 0.75 μm x 2.5 μm, from which Φ = 0.14% can be calculated, and taking into account the sucrose and Ca²⁺ ions in the remediation solution, the total relative viscosity can be expressed as follows, with the final values shown in table 1.

$\mu_r=\mu_{r1}\cdot\mu_{r2}$

Under the constraint of mass conservation, the dynamic behaviour of the fluid satisfies the following conditions:

$\nabla\cdot\vec{u}=0$

Normally, the kinetic behaviour of fluids is analysed by applying numerical simulations, but considering that: 1. this project is more interested in whether the flow process can achieve hypoxia condition; 2. the transformation of energy at the microscopic scale should not only take into account kinetic and potential energy, but also the influence of energy such as the surface energy of the fluid, scale effects will come into play and the situation is more complex; 3. the flow behaviour of the restoration fluid should ideally be consistent with the subsequent The flow behaviour of the repair fluid should ideally be consistent with the subsequent curing behaviour. In summary, we decided to use a cellular automaton approach to simulate the flow of the restoration fluid in the fracture, taking into account both the kinetic behaviour of the fluid and the microscale effects:

To simplify the problem, three particular examples of single fracture models were chosen to consider the problem, and the three single fracture models are shown below:

Fig 2 Three particular examples of single fracture models

In the parallel plate model, the space of the void is set to 4 mm x 15 mm x 2 mm and the side length of each cell is 0.5 mm, the length of the upper convex part in the second model is 5 mm and the height is 1 mm, and the length of the lower convex part in the third model is 5 mm and the height is 0.5 mm. These three simple models allow the distribution of the solid phase in the fracture to be taken into account more efficiently and are therefore used to analysis of hypoxia realisation and flow processe:

In the cellular automaton, the following conditions will be considered and rules developed accordingly:

1. Visible mobility behaviour:

In the flow between plates, the pressure difference of the liquid can be neglected and the remaining external forces are considered only for the effect of gravity and discussed subsequently, the one-dimensional NS equation can be reduced to:

$\frac{\partial u}{\partial t}+\vec{u}\cdot\frac{\partial u}{\partial x}=\frac{\mu}{\rho}\cdot\frac{\partial^2 u}{\partial x^2}$

Given that the behaviour of the restoration fluid flowing from the liquid phase to the void is the case we are mainly considering, and taking into account the discretization of the space by the cellular automaton, this differential equation can be reduced to:

$\frac{du}{dt}=\frac{1}{x}\cdot u^2-\frac{\mu}{\rho x^2}\cdot u(\frac{\partial u}{\partial x}\rightarrow\frac{-u}{x},\frac{\partial^2 u}{\partial x^2}\rightarrow\frac{-u}{x^2})$

x is the edge length of each cell and can be solved for u(t) as

$u=\frac{\mu^2}{\rho\mu x+\rho^2x^4C_0e^{\frac{\mu}{\rho x^2}\cdot t}}$

$C_0=\frac{\mu^2-\rho\mu xu_0}{\rho^2 x^4 u_0}$

u₀ denotes the initial velocity at a given void when there is a fresh inflow of restoration fluid. In a three-dimensional cellular automaton, we consider in this judgement condition only the flow of liquid in a given layer, i.e. a two-dimensional flow problem, and include the isotropic consideration that the liquid flows in both the x and y directions at a velocity u. Thus, for a given cell in the void state, if it is surrounded by water at the moment kT (k ∈ Z), then after the moment T_c1 = [x/u_kT]·T (due to the cellular automaton discretization of time, which will be rounded upwards to ensure that T_c1 is non-zero), the tuple will move from the gap state to the watered state as follows.

$T_{c1}=[\frac{\rho\mu x^2+\rho^2x^5C_0e^{\frac{\mu}{\rho x^2}\cdot kT}}{\mu^2}],k\in Z$

T was chosen to be the same as for the subsequent quorum sensing module, at 0.72s.

2. The effect of gravity on the flow of the repair fluid.

The treatment remains as above, and after taking into account gravity, the NS equation can be reduced to:

$\frac{du}{dt}=\frac{1}{x}\cdot u^2-\frac{\mu}{\rho x^2}\cdot u-g$

To simplify the calculation, we choose a positive equilibrium point of u (du/dt=0) as u_max and linearise it to obtain T_c2=[2x/u_max]·T (rounded upwards). Then, for a given cell in the void state, if there is water above it, the cell changes to the water state and the cell above it to the void state after the moment T_c2.

3. Fluid behaviour at the microscopic scale.

When the velocity of the restoration fluid is low (below a given u_min), the fluid will still behave in the fracture despite not being subject to gravity, and the NS equation still applies at this point, with the term that plays a major role being:

$-\nabla p+\mu\cdot\Delta u=0$

The expression for the pressure difference at this point can be given by the Young-Laplace equation as:

$\Delta p=2H\cdot\gamma$

H is the average curvature of the liquid surface and γ is the surface tension coefficient of the liquid. Although the mean curvature is related to properties such as the contact angle between the solid and liquid phases, the actual calculation is complex and can only give guidance on an order of magnitude basis. We have therefore decided to use a stochastic approach to characterise the process, which is also consistent with the discrete nature of cellular automation.

Gaussian processes are extremely classical stochastic processes, and we consider the width and height of the crack as factors influencing the flow behaviour of the fluid at the microscopic scale, and since the smaller the scale, the stronger the influence of the microscopic conditions, a certain probability function is defined as follows:

$p(x,z)=\frac{1}{\sqrt{2\pi}}\cdot e^{-\frac{1}{2}[\frac{1}{(k_x\mu\rho x)^2}]+\frac{1}{(k_z\mu\rho z)^2}]}$

x denotes the width of a cell, z denotes the height of a cell, and kx and kz are scale-dependent coefficients. Then, for a cell in the void state at a certain location, the probability of its transition to the watery state after 5T is related to the number of surrounding four neighbours in the watery state, N, and can be calculated as follows.

$P_{empty\rightarrow water}=1-(1-\int\int p(x,z)dxdz)^N$

The hypoxia case is considered as a proportion of the remaining void space to the void space at the initial moment and is expressed as follows.

$vocuum=\frac{V(kT)}{V_{initial}}$

When the vocuum is < 0.1, it is considered to be in a hypoxia state. The parameters used in this section are listed in the following table.

Table 1 Parameters in flow part

The flow of the restoration fluid can be divided into two parts in this cellular automation model, initially the macroscopic flow behaviour, when [x/u] > 100T, at which point the stochastic flow under microscale effects is considered and the simulation ends after a further 400T.

The flow simulation processes for the three single fracture models under the above conditions and the variation in void fraction are shown in Figure 3:

Fig 3 (a) Flow in the parallel plate fracture, entering the hypoxia state after 340T (b) Flow in the upper convex model, entering the hypoxia state after 330T (c) Flow in the lower convex model, entering the hypoxia state after 330T

All three simple single fracture models show that the restoration fluid can be brought to a hypoxia state after a certain period of time and that the restoration fluid can occupy the entire crack well after a longer period of time. Therefore, in more complex fractures, a hypoxia state can also be reached in most cases to provide a basis for the initiation of subsequent modules, and the flow conditions will also be well prepared for the solidification process. In the subsequent modelling section, only the parallel plate model is used as the basis for the model.

Also, to account for the concentration of the substance in the configuration of the restoration fluid, we set the initial viscosity and density to a series of values with gradient differences and expressed the void ratio for each case at 300T and 400T after the start of random flow is shown in Figure 4:

Fig 4 Proportion of voids at 300T and 400T for different initial conditions

From the bar graph we can conclude that hypoxia conditions can be achieved at the end of the stochastic flow process when the viscosity and density are small, and the previous equation gives the viscosity and density versus added substances and bacteria, so we can give control of substance addition based on the constraints of low oxygen conditions.

The theoretical basis for the model is derived mainly from [5], in which the normalised population level X₁, trophic level X₂, X₃=[SigH], X₄=[SpaRK], X₅=[SpaS] form a continuum X=[X₁,X₂,X₃,X₄,X₅]^T that together describe the state of the system. s=[S₃,S₄,S₅]^T consists of three binary switches S=[S₃,S₄,S₅]^T consists of three binary switches representing whether the three genes SigH, SpaRK, and SpaS are in the activated (=1) or repressed (=0) state, respectively.

The change of state of the continuous part can be represented by a deterministic model as follows:

$\dot{x_1}=rx_1(1-\frac{x_1}{D_{\infty}(x_2)})$

$\dot{x_2}=-k_1x_1+k_2x_5$

$x_i(t)=\left \{ \begin{array}{**lr**}-l_ix_i(t), & S_i(t)=0, &\\ -l_ix_i(t)+k_i, & S_i(t)=1,&\end{array} \right. i=3,4,5$

The Bacillus subtilis population can be represented by a logistic model where changes in trophic levels are regulated by population and subtilin levels (equating SpaS with mature subtilin), while changes in SigH, SpaRK and SpaS are influenced by the discrete state Si, where the activation state of the genes determines whether they will be produced outside natural degradation. In models of population change, environmental accommodation is a parameter determined by trophic level, defined here as:

$D_{\infty}(x_2)=min\{\frac{x_2}{X_0},D_{max}\}$

X₀, D_max are given as constants. The description of the SigH production state is determined by the nutrient level and is given below:

$S_i(t)=\left\{ \begin{array}{**lr**}1,& x_2(t)\lt n&\\0,& x_2(t) \geq n &\end{array} \right.$

The threshold n indicates that the chymotrypsin production mechanism is excited below this condition. In [6], Markov chains based on probability transfer matrices, the selection of coefficients follows inherent biochemical constraints.

$P[S_i(kT+T)=1|S_i(kT)=0]=\frac{c_i x_{i-1}(kT)}{1+c_i x_{i-1}(kT)}$

$P[S_i(kT+T)=0|S_i(kT)=1]=\frac{1}{1+c_i x_{i-1}(kT)}$

In our design, another major trigger for excitation of chymotrypsin production is hypoxia conditions, so we consider introducing the discrete state of whether hypoxia conditions excite as S6 in the original model. The excitation of low oxygen conditions can be described by the vocuum < 0.1 given in the flow model, expressed as follows.

$S_6(t)=\left \{ \begin{array}{**lr**}1,& vocuum\lt 0.1 & \\0,& vocuum\geq0.1 & \end{array}\right.$

After the hypoxia condition is reached, SigH is skipped and the initiation of SpaRK is directly stimulated. To achieve the effect of this hypoxia condition and to introduce a certain degree of randomness, S4 and S5 are given the following form of change after the introduction of this condition.

$S_i(kT)=1, rand>(1-vocuum)\times(S_6(kT)==1),i=4,5$

The method of discretizing deterministic processes is given in [5], for the macroscopic level x1, x2:

$x_1(kT+T)=\frac{D_{\infty}(x_2)}{1+(\frac{D_{\infty}(x_2)}{x_1}-1)\tilde{r}}$

$x_2(kT+T)=x_2(kT)-\tilde{k_1}x_1(kT)+\tilde{k_2}x_5(kT)$

For the micro level x3, x4, x5:

$x_i(kT+T)= \left \{ \begin{array}{**lr**}\tilde{l_i}x_i(kT),& S_i(kT)=0 & \\ \tilde{l_i}x_i{kT}+\tilde{k_i},& S_i(kT)=1 \end{array}\right.i=3,4,5$

The coefficients of the discretization process are expressed as follows:

$\tilde{r}=e^{-rT}$

$\tilde{k_1}=Tk_1$

$\tilde{k_2}=Tk_2$

$\tilde{l_i}=e^{-l_iT}(i=3,4,5)$

$\tilde{k_i}=k_i\frac{1-\tilde{l_i}}{l_i}(i=3,4,5)$

The parameters used in the model are listed below:

Table 2 parameters in quorum sensing module

The variation of the states of five substances with time is shown in Figure 5:

Fig 5 Normalisation of population, nutrient, sigH, spaRK and spaS in the hybrid model. and spaS, with the moments when spaS is activated to a certain level

The moment of spaS promoter initiation was more influenced by hypoxia conditions in the simulation. After excitation by hypoxia conditions, the spaS promoter intensity reached half of its peak at 359T (in the experiment, 0.1 mg/ml of spaS was obtained to trigger downstream expression, and it was considered that half of the peak had been met), at which point the CA module and ACCBP module were considered to have started.

The Carbonic Anhydrase Production is responsible for the output of carbonate ions. Carbonic anhydrase can catalyse the conversion of CO₂ to HCO₃^-in the presence of Zn²⁺. The HCO₃^- that is transported out of the cell down the concentration gradient is combined with Ca²⁺ under alkaline conditions to produce amorphous calcium carbonate, the biochemical process of which can be expressed as follows:

$DNA\rightarrow mRNA_{CA}$

$mRNA_{CA}\rightarrow CA$

$CO_2+H_2O\xrightarrow{CA/Zn^{2+}}HCO_3^-$

$HCO_3^-+Ca^{2+}+OH^-\rightarrow CaCO_3+H_2O$

Both the production and degradation of intracellular macromolecules are considered as first-order linear kinetic reactions, which can be expressed as follows:

$\frac{d[mRNA_{CA}]}{dt}=\alpha_1\cdot[DNA]-\beta_1\cdot[mRNA_{CA}]$

$\frac{d[CA]}{dt}=\alpha_2\cdot[mRNA_{CA}]-\beta_2\cdot[CA]$

[DNA] denotes the number of plasmid copies in each Bacillus subtilis, estimated to be 10, α₁ denotes the transcription rate of the gene and α₂ integrates the processes of translation of the mRNA, folding of the polypeptide, etc. The average elongation rates of transcription and translation in B. subtilis are given by [7] as 73nt/s and 47nt/s, respectively, and it was observed that most operons in B. subtilis are transcribed without a closely trailing ribosome, and our gene length is 837 bp, from which α₁ can be calculated. Also the total elapsed time for the protein folding process can be obtained from [8], and combined with the translation elapsed time, α₂ can be estimated.

The half-life of CA degradation is given in [9] as 8 h. The half-life of mRNA in Bacillus subtilis is illustrated in [10] as 30s~more than 20 minates, and conservatively, the half-life of mRNA_CA was selected as 10 min. According to the first-order linear equation of the degradation process, we can obtain:

$\frac{dx}{dt}=-\beta x\rightarrow T_{\frac{1}{2}}=\frac{ln2}{\beta}$

From this, β₁ and β₂ can be estimated. The chemical reactions involved in carbonic anhydrase have been studied in [11], and in the carboxylation reaction of carbon dioxide, taking into account the contribution of both catalytic and non-catalytic components, the rates of each reaction can be expressed as follows (assuming that all of the HCO₃^- produced can be transported outside the cell):

$\frac{dc[HCO_3^-]}{dt}=(k_{+2}+\frac{k_{Cat}}{K_M}\cdot c[CA])\cdot c[CO_2]$

$\frac{dc[CaCO_3]}{dt}=K_S\cdot c[HCO_3^-]\cdot c[OH^-]\cdot c[Ca^{2+}]$

To illustrate that the Carbonic Anhydrase Production we have used serves to accelerate the process of calcium carbonate precipitation, the determination of the rate-limiting step needs to be discussed. The rate of carbon dioxide production in Bacillus subtilis, k_CO2, can be traced in the graph of [12] as 5-20 mmol·g^-1·h^-1, and its consumption rate is the rate of HCO₃^- production. Thus, without the use of CA enzyme catalysis, the rate of CO₂ production and consumption can be expressed as:

$v_{CO_2+}=k_{CO_2}$

$v_{CO_2-}=k_{+2}\cdot c_0[CO_2]$

Where c₀[CO₂] is the concentration of carbon dioxide in water (i.e. the solubility of carbon dioxide in water), it can be calculated that v_CO2+ > v_CO2-, so that the bicarbonate production step in the precipitation of calcium carbonate is the rate-limiting step and the CA enzyme accelerates the process very well.

In chemical reactions, each reactant unit is a concentration, so the [CA] quantity in a biomolecular reaction needs to be converted to a concentration in the following form:

$c[CA]=\frac{[CA]\cdot c_{Bacillus}}{N_A}$

c_N denotes the concentration of Bacillus subtilis in the restoration solution, which is 10⁸ per mL and is given in the flow procedure; c[CA] is the fraction of the substance of carbonic anhydrase in the remediation solution.In summary, the reaction rates involved in the chemical reaction section can be given as follows:

$\frac{dc[HCO_3^-]}{dt}=(k_{+2}+\frac{k_{Cat}}{K_M}\cdot\frac{c_{Bacillus}}{N_A}c[CA])\cdot c[CO_2]$

$\frac{dc[CaCO_3]}{dt}=K_S\cdot c[HCO_3^-]\cdot c[OH^-]\cdot c[Ca^{2+}]$

$\frac{dc[CO_2]}{dt}=k_{CO_2}-(k_{+2}+\frac{k_{Cat}}{K_M}\cdot\frac{c_{Bacillus}}{N_A}[CA])\cdot c[CO_2]$

$\frac{dc[Ca^{2+}]}{dt}=-K_S\cdot c[HCO_3^-]\cdot c[OH^-]\cdot c[Ca^{2+}]$

The parameters involved in the module and their significance are indicated in the following table:

Table 3 parameters in Carbonic Anhydrase Production

The variation of mRNA_CA, [CA], c[HCO₃^-], c[CaCO₃] with time is shown in Figure 6 with a statistical duration of 36000s, i.e. 10h.

Fig 6 Changes in mRNA_CA, [CA], c[HCO₃^-], c[CaCO₃] with time

The end product CaCO₃ over time will be saved and used in the solidification model.

The Biogolical Scaffold Module is initiated at the same time as the Carbonic Anhydrase Production. The naturally precipitated calcium carbonate is loosely structured, whereas the biosynthetic bioscaffold has a good anchoring effect on the amorphous calcium carbonate. The Biogolical Scaffold Module can produce proteins such as spytag, spycatcher and EutM to combine with the bacterial cell wall to form the bioscaffold, and the biochemical process can be described as follows:

$DNA\rightarrow mRNA_{EutM-Spycatcher}$

$DNA\rightarrow mRNA_{spytag-Hag}$

$mRNA_{EutM-Spycatcher}\rightarrow EutM-Spycatcher$

$mRNA_{spytag-Hag}\rightarrow spytag-Hag$

$EutM-Spycatcher\rightarrow (EutM-Spycatcher)_n-complex1$

$complex1+spytag-Hag\rightarrow complex2$

The deterministic process can be described by a system of differential equations as follows:

$\frac{d[mRNA_{EutM-Spycatcher}]}{dt}=\alpha_3\cdot DNA-\beta_1\cdot [mRNA_{EutM-Spycatcher}]$

$\frac{d[mRNA_{spytag-Hag}]}{dt}=\alpha_4\cdot DNA-\beta_1\cdot [mRNA_{spytag-Hag}]$

$\frac{d[EutM-Spycatcher]}{dt}=\alpha_5\cdot[mRNA_{EutM-Spycatcher}]-\beta_3\cdot[EutM-Spycatcher]-\gamma_1\cdot[EutM-Spycatcher]+\mu_1\cdot[complex1]$

$\frac{d[spytag-Hag]}{dt}=\alpha_6\cdot[mRNA_{spytag-Hag}]-\beta_4\cdot[spytag-Hag]-\gamma_2\cdot[complex1]\cdot[spytag-Hag]+\mu_1[complex2]$

$\frac{d[complex1]}{dt}=\gamma_1\cdot[EutM-Spycatcher]-\mu_1\cdot[complex1]+\mu_2\cdot[complex2]-\beta_5\cdot[complex1]$

$\frac{d[complex2]}{dt}=\gamma_2\cdot[spytag-Hag]\cdot[complex1]-\mu_2\cdot[complex2]-\beta_6\cdot[complex2]$

α₃ and α₄ indicate the transcription rate of the gene. The gene lengths of EutM-Spycatcher and Spytag-Hag are 771bp and 954bp respectively, and their transcription rates can be calculated by taking the same approach as in the Quorum Sensing Module; α₅ and α₆ indicate the rate of translation, which is not found in the literature and takes the value of the order of magnitude in the CA module; β₁ is the rate of mRNA degradation, the same as in the CA module; β₃, β₄, β₅, β₆ are the rates of degradation of each complex protein, the half-life of FNBP is given by [13] as 46 min, and its degradation rate can be obtained by referring to this; γ₁, γ₂ denote the rate of protein binding, the value of γ₂ is given by [14], μ₃, μ₄ are the rates of dissociation of the complex. The values of γ₁, μ₃ and μ₄ were not obtained from the relevant literature and were chosen as empirical values for the time being, the effects of which will be discussed in the subsequent sensitivity analysis.

The parameters involved in this part and their significance are indicated in the following table:

Table 4 parameters in Biogolical Scaffold Module

The non-rigid ODE solver is used to solve the above system of equations and the end product complex2, i.e. the change in the number of biological scaffolds, is represented in Figure 7 with a statistical time of 36,000s, i.e. 10h.

Fig 7 the change of scaffolds over time

Changes in Scaffolds over time will be saved and eventually used in the curing process.

The reliability of the selection of values for γ₁, μ₃ and μ₄ can be analysed in terms of the sensitivity of these three parameters, by comparing their values with the original results by scaling them up and down by a factor of ten respectively and plotting them in Figure 8.

Fig 8 (a) Effect of γ₁ on scaffold concentration at ten times magnification and reduction respectively (b) Effect of μ₃ and μ₄ on scaffold concentration at ten times magnification and reduction respectively

From Figure 8 it can be seen that γ₁, μ₃ and μ₄ have little or no effect on the variation in the number of final scaffolds, so the empirical values chosen for these three parameters are plausible.

The initial state of the solidificaiton process uses the final state of the parallel plate model flow used in the flow process, in which the simulation of the solidification of the repair fluid is carried out by means of certain rules.

The work in [16] shows that the solidification behaviour under the MICP method starts at the surface, and we will start with this feature to give the corresponding curing rules.

1. Solidification substance concentration.

In our design, the main solidification material is amorphous calcium carbonate produced by the Carbonic Anhydrase Production, in addition to a composite of spytag-spycatcher produced by the ACCBP module, EutM, bacterial cell walls and other scaffolds.

In the carbonate anhydrase module and Biogolical Scaffold Module the changes in CaCO₃ and Scaffolds have been obtained over 36000 s. During the curing process the time interval of change T = 360 s is set. Each T moment the concentration of CaCO₃ and Scaffolds concentration was chosen to be the value that exceeded that moment and was the closest.

Therefore, we believe that if the two main solidification substances, amorphous calcium carbonate and biological scaffold, reach a certain concentration of n_ca and n_scaffold in a certain aqueous state of the cell near the surface of the solid, the cell will be converted to a cured state of the cell with a certain probability P_{solidifacation} at the next moment, which is expressed as follows.

$P_{Solidification}=\frac{c[CaCO_3(kT)]\cdot[scaffolds](kT)}{c[CaCO_3(kT)]_{max}\cdot[scaffolds]_{max}},c[CaCO_3(kT)]\gt 0.0025,[scaffolds]\gt 2000$

2. Destruction of new curing material

We consider that the newly generated solidification is not necessarily stable in the presence of surrounding restoration fluids, so that for a given solidified cell, the greater the number of surrounding neighbours in the aqueous state, the easier it is for the cell to be reconverted to the aqueous state, and the probability of this process occurring is expressed as follows.

$P_{solid\rightarrow water}=1-(1-p_{Solidification}^N)$

N denotes the total number of water-bearing neighbours around the cell.

The results of the simulation of the solidification process can be represented as follows (Figure 9), where the yellow colour is the newly created solidification material.

Fig 9 The simulation of the solidificaiton process

It can be seen that the final solidification result is more than satisfactory.