Diffusion/Reaction Limited Aggregation Approach for Microstructure Evolution and Condensation Kinetics during Synthesis of Silica-Based Alcogels

A base-catalysed methyltrimethoxysilane (MTMS) colloidal gel formation was implemented as a cellular automaton (CA) system, specifically diffusion and/or reaction-limited aggregation. The initial characteristic model parameters were determined based on experimental synthesis of MTMS-based, ambient-pressure-dried aerogels. The applicability of the numerical approach to the prediction of gels’ condensation kinetics and their structure was evaluated. The developed model reflects the kinetics properly within the investigated chemical composition range (in strongly reaction-limited aggregation conditions) and, to a slightly lesser extent, the structural properties of aggregates. Ultimately, a relatively simple numerical model reflecting silica-based gel formation was obtained and verified experimentally. The CA simulations have proved valid for understanding the relation between the initial chemical composition and kinetics constants of MTMS-based synthesis and their impact on secondary particle aggregation process kinetics.


Introduction
The sol-gel synthesis method is based on a transition between two colloidal systems. The transition may proceed due to hydrolysis, poly-condensation and re-esterification reactions [1][2][3]. The reactions occur, to a certain pH-dependent extent, concurrently [4,5]. It is a versatile method of producing many industrially and commercially important materials such as ceramics [6], glasses [7], coating films [8,9], adsorbents [10], insulation materials [11] and catalyst supports [3,12,13]. The chemistry of the transition strongly depends on the chosen precursors and the initial chemical composition of a reaction mixture. The final product's structure-film, xerogel, aerogel, ceramic fibres, dense ceramics-depends on the chosen preparation and drying method [13].
During gel formation from organoalkoxysilane precursor molecules, three fundamental steps can be distinguished-the first two have already been mentioned: (i) hydrolysis and (ii) condensation, while (iii) is phase separation [2]. Phase separation was described by Issa and Luyt as a process when a reaction mixture loses its homogeneity to form a liquid in a continuous solid phase or a solid in a continuous liquid phase [2]. The first case is the definition of a gel, while the second is a sol or a precipitate (in the case of sedimentation). Mechanisms of phase separation for silica precursors were exhaustively studied by Nakanishi et al. [14][15][16]. The imposition of different chemical reactions and physical phenomena makes the investigation of organoalkoxysilanes formation kinetics a challenging yet necessary field of study. With such a variety of simultaneous factors that are impossible to investigate separately experimentally, the combination of modelling and experimental studies could provide an in-depth understanding of the fundamentals of the sol-gel transition. Studies on gel formation kinetics are also an essential basis for process optimisation. and experimental studies could provide an in-depth understanding of the fundamentals of the sol-gel transition. Studies on gel formation kinetics are also an essential basis for process optimisation.
Every method is associated with different limitations. The Smoluchowski equation does not provide structural information and is limited to a kinetics description [43]. Brownian and Langevin dynamics require a knowledge/designation of peculiar physicochemical parameters (based on our experience with coarse-grained molecular dynamics studies of polyelectrolytes' pH-dependent conformations [44]). Finally, percolation theory and DLA/RLA models were initially recognised as approaches that provided a good view of the internal structure but no direct information about system evolution kinetics [31,41]. However, this disadvantage of DLA/RLA was overcome, e.g., by Hsieh et al. [41] and in our previous work [41,45].
These approaches can be implemented via different models of computation, depending on the chosen scale ( Figure 1), such as the Monte Carlo method [46], molecular dynamics (reactive, coarse-grained) [1,[36][37][38][47][48][49], population balance [50] and cellular automata. The last one is rarely applied for gelation description [42], though is more common for investigating mechanical behaviour [51,52]. However, cellular automata are used for many processes concerning aggregation due to diffusion (possibly with chemical reaction) [53], such as recrystallisation [54] or flocculation [55]. The most significant challenges associated with modelling studies are compliance with experimental observations, model utility and scope of applicability. In the presented paper, we focused on providing an approach to the translation of the synthesis parameters of methyltrimethoxysilane (MTMS) gels on characteristic parameters in a developed 3D DLCA/RLCA cellular automaton. CA is a well-known method, although most studies lack any thoughts on particle/cluster connection probability and only concern diffusion-limited aggregation [41,42]. The presented work focuses on the interpretation of this The most significant challenges associated with modelling studies are compliance with experimental observations, model utility and scope of applicability. In the presented paper, we focused on providing an approach to the translation of the synthesis parameters of methyltrimethoxysilane (MTMS) gels on characteristic parameters in a developed 3D DLCA/RLCA cellular automaton. CA is a well-known method, although most studies lack any thoughts on particle/cluster connection probability and only concern diffusion-limited aggregation [41,42]. The presented work focuses on the interpretation of this parameter and finding the relation between experimental determinable values, which is the main novelty advancing the state of the art.

Estimation of Model Parameters
Activation energy E and A constant were established based on the conducted measurements of gelation time for A-F syntheses for the chosen temperatures (described in Section 4.1). The results were plotted in ln t g = f (1/T) form and presented in Figure 2, along with the fitted linear trendlines (R 2 values are presented in the graph). According to the linearisation of the Arrhenius equation (Equations (1) and (2)), the trendline constants allowed for an estimation of the E and A values for the investigated syntheses (presented in Figure 3). ln A denotes the Arrhenius constant, R-gas constant, and E-activation energy. A is a constant dependent on the Arrhenius constant. The A constant can be interpreted as the number of molecular collisions with proper orientation, and the exponential element exp − E RT g as the probability of a successful collision; therefore, k is the number of effective molecular collisions per second.
denotes the Arrhenius constant, -gas constant, and -activation energy.
is a constant dependent on the Arrhenius constant. The constant can be interpreted as the number of molecular collisions with proper orientation, and the exponential element exp − as the probability of a successful collision; therefore, is the number of effective molecular collisions per second. These values were a basis for the probability parameter designation (according to the interpretation of the Arrhenius equation), which is an equivalent of an effective collision probability value and the primary model parameter within the presented research. The probability rate was designated directly based on the Arrhenius equation, while probability was multiplied by the time step value (characteristic for each synthesis) (Equation (4)). The probability and probability rate dependence on the catalyst/precursor molar ratio are shown in Figure 4.
Timestep for each case was calculated on the basis of our previous research; we have experimentally measured (using the tilting test-tube method) the value of gelation time, and we can interpolate a predicted numerical value of a gel-point based on the interpolation of 25 studied cases.    These values were a basis for the probability parameter designation (according to the interpretation of the Arrhenius equation), which is an equivalent of an effective collision probability value and the primary model parameter within the presented research. The probability rate was designated directly based on the Arrhenius equation, while probability P was multiplied by the time step value (characteristic for each synthesis) (Equation (4)). The probability and probability rate dependence on the catalyst/precursor molar ratio are shown in Figure 4.
Timestep for each case was calculated on the basis of our previous research; we have experimentally measured (using the tilting test-tube method) the value of gelation time, and we can interpolate a predicted numerical value of a gel-point based on the interpolation of 25 studied cases.  The cellular automaton requires two set parameters: (i) secondary particle concentration and (ii) probability of an effective collision.
Within the numerical studies, a set of 25 cases was simulated: = 0.5, 1, 2, 3, 4 [% ] and = 10 , 10 , 10 , 10 , 1 [−] . For example, a concentration of 1% means that secondary particles occupy 1% of cells. As to the probability value, = 1 means there is a 100% probability that two adjacent particles (located in adjacent cells, according to Moore's neighbourhood) will become an aggregate (thus, the DLA mechanism is responsible for the aggregation). For ≪ 1, most collisions do not end in a permanent association of particles/aggregates (RLA mechanism). Subsequently, the experimental conditions had to be "translated" into model parameters. The first step was an estimation of the secondary particle concentration in a sample. Knowing the volume of an alcogel, the mass of the skeleton (aerogel mass), and structural density, the volume concentration of secondary particles in a system was calculated according to Equation (1): where is the volume of secondary particles, is the volume of alcogel, is thenumber of secondary particles in the system, is the volume of one particle, m_s is the mass of skeleton (aerogel mass), and is structural density (1.38 [g cm ⁄ ]). The cellular automaton requires two set parameters: (i) secondary particle concentration and (ii) probability of an effective collision.
Within the numerical studies, a set of 25 cases was simulated: c = {0.5, 1, 2, 3, 4 [% vol ]} and P = 10 −4 , 10 −3 , 10 −2 , 10 −1 , 1 [−] . For example, a concentration of 1% means that secondary particles occupy 1% of cells. As to the probability value, P = 1 means there is a 100% probability that two adjacent particles (located in adjacent cells, according to Moore's neighbourhood) will become an aggregate (thus, the DLA mechanism is responsible for the aggregation). For P 1, most collisions do not end in a permanent association of particles/aggregates (RLA mechanism). Subsequently, the experimental conditions had to be "translated" into model parameters. The first step was an estimation of the secondary particle concentration in a sample. Knowing the volume of an alcogel, the mass of the skeleton (aerogel mass), and structural density, the volume concentration of secondary particles in a system was calculated according to Equation (1): where V sp is the volume of secondary particles, V a is the volume of alcogel, n sp is thenumber of secondary particles in the system, V sp1 is the volume of one particle, m s is the mass of skeleton (aerogel mass), and ρ s is structural density (1.38 g/cm 3 ). Probability was evaluated based on the Arrhenius equation, but recounted per timestep dt, because this value is the most adequate to the model parameter (Equation (3)).
Timestep for each case was calculated on the basis of our previous research; we have experimentally measured (using the tilting test-tube method) the value of gelation time, and we can interpolate a predicted numerical value of a gel-point based on the interpolation of 25 studied cases.
All of the designated model parameters, c, P and dt for actual experimental syntheses A-F, are gathered in Table 1. Table 1. Designated model parameters: secondary particle concentration (c), probability of an effective collision (P) and time step (dt) for actual experimental syntheses A-F. timestep dt, because this value is the most adequate to the model parameter.
Timestep for each case was calculated on the basis of our previous research; we have experimentally measured (using the tilting test-tube method) the value of gelation time, and we can interpolate a predicted numerical value of a gel-point based on the interpolation of 25 studied cases.
All of the designated model parameters, c, P and dt for actual experimental syntheses A-F, are gathered in Table 1. Table 1. Designated model parameters: secondary particle concentration (c), probability of an effective collision (P) and time step (dt) for actual experimental syntheses A-F.

Kinetics Analysis
While analysing the kinetics curves, three phases can be distinguished. The first one refers to nucleation, which happens on the molecular level, and does not cause the turbidity increase. The solution is clear, and mass growth is negligible. The second phase is visible as an intensive increase in absorbance value that can be correlated with silica skeleton formation via siloxane bonding. A plateau informs us of reaching the gel point where gel is in dynamic equilibrium-gel can be reorganised via coarsening and Ostwald ripening.
Within the model, the first phase (nucleation) cannot be observed in our numerical system-particles in CA are already formed, and their aggregation leads directly to the second phase. The third phase, regarded as a plateau, is when mostly aggregate clusters that were already detectable impact the structure formation process and morphology, but does not translate into an increase in the mass of the condensing product.
A summary of these theoretical considerations and assumptions is gathered in Table  2. We identified condensation mechanisms in each of the distinguished phases, according to Figure 5. During the first phase, the mass of the condensing product does not (or barely) increase. It means that the hydrolysed precursor molecules' condensation leads to nuclei formation, but they are too small to be detected by UV-Vis spectrophotometer. During the second phase, we have intensive mass growth. It means that large aggregates are becoming rapidly detectable. A mechanism that would explain such phenomena is bonding two undetectable aggregates into a bigger one which is immediately registered. Another option (however, less rapid) is when a small undetectable aggregate bonds with c = 6.7% P = 1.4 × 10 −4 dt = 1.34 × 10 −1 timestep dt, because this value is the most adequate to the model parameter.
Timestep for each case was calculated on the basis of our previous research; we have experimentally measured (using the tilting test-tube method) the value of gelation time, and we can interpolate a predicted numerical value of a gel-point based on the interpolation of 25 studied cases.
All of the designated model parameters, c, P and dt for actual experimental syntheses A-F, are gathered in Table 1. Table 1. Designated model parameters: secondary particle concentration (c), probability of an effective collision (P) and time step (dt) for actual experimental syntheses A-F.

Kinetics Analysis
While analysing the kinetics curves, three phases can be distinguished. The first one refers to nucleation, which happens on the molecular level, and does not cause the turbidity increase. The solution is clear, and mass growth is negligible. The second phase is visible as an intensive increase in absorbance value that can be correlated with silica skeleton formation via siloxane bonding. A plateau informs us of reaching the gel point where gel is in dynamic equilibrium-gel can be reorganised via coarsening and Ostwald ripening.
Within the model, the first phase (nucleation) cannot be observed in our numerical system-particles in CA are already formed, and their aggregation leads directly to the second phase. The third phase, regarded as a plateau, is when mostly aggregate clusters that were already detectable impact the structure formation process and morphology, but does not translate into an increase in the mass of the condensing product.
A summary of these theoretical considerations and assumptions is gathered in Table  2. We identified condensation mechanisms in each of the distinguished phases, according to Figure 5. During the first phase, the mass of the condensing product does not (or barely) increase. It means that the hydrolysed precursor molecules' condensation leads to nuclei formation, but they are too small to be detected by UV-Vis spectrophotometer. During the second phase, we have intensive mass growth. It means that large aggregates are becoming rapidly detectable. A mechanism that would explain such phenomena is bonding two undetectable aggregates into a bigger one which is immediately registered. Another option (however, less rapid) is when a small undetectable aggregate bonds with c = 6.4% P = 1.8 × 10 −4 dt = 1.31 × 10 −1 timestep dt, because this value is the most adequate to the model parameter.
Timestep for each case was calculated on the basis of our previous research; we have experimentally measured (using the tilting test-tube method) the value of gelation time, and we can interpolate a predicted numerical value of a gel-point based on the interpolation of 25 studied cases.
All of the designated model parameters, c, P and dt for actual experimental syntheses A-F, are gathered in Table 1. Table 1. Designated model parameters: secondary particle concentration (c), probability of an effective collision (P) and time step (dt) for actual experimental syntheses A-F.

Kinetics Analysis
While analysing the kinetics curves, three phases can be distinguished. The first one refers to nucleation, which happens on the molecular level, and does not cause the turbidity increase. The solution is clear, and mass growth is negligible. The second phase is visible as an intensive increase in absorbance value that can be correlated with silica skeleton formation via siloxane bonding. A plateau informs us of reaching the gel point where gel is in dynamic equilibrium-gel can be reorganised via coarsening and Ostwald ripening.
Within the model, the first phase (nucleation) cannot be observed in our numerical system-particles in CA are already formed, and their aggregation leads directly to the second phase. The third phase, regarded as a plateau, is when mostly aggregate clusters that were already detectable impact the structure formation process and morphology, but does not translate into an increase in the mass of the condensing product.
A summary of these theoretical considerations and assumptions is gathered in Table  2. We identified condensation mechanisms in each of the distinguished phases, according to Figure 5. During the first phase, the mass of the condensing product does not (or barely) increase. It means that the hydrolysed precursor molecules' condensation leads to nuclei formation, but they are too small to be detected by UV-Vis spectrophotometer. During the second phase, we have intensive mass growth. It means that large aggregates are becoming rapidly detectable. A mechanism that would explain such phenomena is bonding two undetectable aggregates into a bigger one which is immediately registered. Another option (however, less rapid) is when a small undetectable aggregate bonds with timestep dt, because this value is the most adequate to the model parameter.
Timestep for each case was calculated on the basis of our previous research; we have experimentally measured (using the tilting test-tube method) the value of gelation time, and we can interpolate a predicted numerical value of a gel-point based on the interpolation of 25 studied cases.
All of the designated model parameters, c, P and dt for actual experimental syntheses A-F, are gathered in Table 1. Table 1. Designated model parameters: secondary particle concentration (c), probability of an effective collision (P) and time step (dt) for actual experimental syntheses A-F.

Kinetics Analysis
While analysing the kinetics curves, three phases can be distinguished. The first one refers to nucleation, which happens on the molecular level, and does not cause the turbidity increase. The solution is clear, and mass growth is negligible. The second phase is visible as an intensive increase in absorbance value that can be correlated with silica skeleton formation via siloxane bonding. A plateau informs us of reaching the gel point where gel is in dynamic equilibrium-gel can be reorganised via coarsening and Ostwald ripening.
Within the model, the first phase (nucleation) cannot be observed in our numerical system-particles in CA are already formed, and their aggregation leads directly to the second phase. The third phase, regarded as a plateau, is when mostly aggregate clusters that were already detectable impact the structure formation process and morphology, but does not translate into an increase in the mass of the condensing product.
A summary of these theoretical considerations and assumptions is gathered in Table  2. We identified condensation mechanisms in each of the distinguished phases, according to Figure 5. During the first phase, the mass of the condensing product does not (or barely) increase. It means that the hydrolysed precursor molecules' condensation leads to nuclei formation, but they are too small to be detected by UV-Vis spectrophotometer. During the second phase, we have intensive mass growth. It means that large aggregates are becoming rapidly detectable. A mechanism that would explain such phenomena is bonding two undetectable aggregates into a bigger one which is immediately registered. Another option (however, less rapid) is when a small undetectable aggregate bonds with timestep dt, because this value is the most adequate to the model parameter.
Timestep for each case was calculated on the basis of our previous research; we have experimentally measured (using the tilting test-tube method) the value of gelation time, and we can interpolate a predicted numerical value of a gel-point based on the interpolation of 25 studied cases.
All of the designated model parameters, c, P and dt for actual experimental syntheses A-F, are gathered in Table 1. Table 1. Designated model parameters: secondary particle concentration (c), probability of an effective collision (P) and time step (dt) for actual experimental syntheses A-F.

Kinetics Analysis
While analysing the kinetics curves, three phases can be distinguished. The first one refers to nucleation, which happens on the molecular level, and does not cause the turbidity increase. The solution is clear, and mass growth is negligible. The second phase is visible as an intensive increase in absorbance value that can be correlated with silica skeleton formation via siloxane bonding. A plateau informs us of reaching the gel point where gel is in dynamic equilibrium-gel can be reorganised via coarsening and Ostwald ripening.
Within the model, the first phase (nucleation) cannot be observed in our numerical system-particles in CA are already formed, and their aggregation leads directly to the second phase. The third phase, regarded as a plateau, is when mostly aggregate clusters that were already detectable impact the structure formation process and morphology, but does not translate into an increase in the mass of the condensing product.
A summary of these theoretical considerations and assumptions is gathered in Table  2. We identified condensation mechanisms in each of the distinguished phases, according to Figure 5. During the first phase, the mass of the condensing product does not (or barely) increase. It means that the hydrolysed precursor molecules' condensation leads to nuclei formation, but they are too small to be detected by UV-Vis spectrophotometer. During the second phase, we have intensive mass growth. It means that large aggregates are becoming rapidly detectable. A mechanism that would explain such phenomena is bonding two undetectable aggregates into a bigger one which is immediately registered. Another option (however, less rapid) is when a small undetectable aggregate bonds with c = 6.1% P = 3.5 × 10 −4 dt = 1.12 × 10 −1 timestep dt, because this value is the most adequate to the model parameter.
Timestep for each case was calculated on the basis of our previous research; we have experimentally measured (using the tilting test-tube method) the value of gelation time, and we can interpolate a predicted numerical value of a gel-point based on the interpolation of 25 studied cases.
All of the designated model parameters, c, P and dt for actual experimental syntheses A-F, are gathered in Table 1. Table 1. Designated model parameters: secondary particle concentration (c), probability of an effective collision (P) and time step (dt) for actual experimental syntheses A-F.

Kinetics Analysis
While analysing the kinetics curves, three phases can be distinguished. The first one refers to nucleation, which happens on the molecular level, and does not cause the turbidity increase. The solution is clear, and mass growth is negligible. The second phase is visible as an intensive increase in absorbance value that can be correlated with silica skeleton formation via siloxane bonding. A plateau informs us of reaching the gel point where gel is in dynamic equilibrium-gel can be reorganised via coarsening and Ostwald ripening.
Within the model, the first phase (nucleation) cannot be observed in our numerical system-particles in CA are already formed, and their aggregation leads directly to the second phase. The third phase, regarded as a plateau, is when mostly aggregate clusters that were already detectable impact the structure formation process and morphology, but does not translate into an increase in the mass of the condensing product.
A summary of these theoretical considerations and assumptions is gathered in Table  2. We identified condensation mechanisms in each of the distinguished phases, according to Figure 5. During the first phase, the mass of the condensing product does not (or barely) increase. It means that the hydrolysed precursor molecules' condensation leads to nuclei formation, but they are too small to be detected by UV-Vis spectrophotometer. During the second phase, we have intensive mass growth. It means that large aggregates are becoming rapidly detectable. A mechanism that would explain such phenomena is bonding two undetectable aggregates into a bigger one which is immediately registered. Another option (however, less rapid) is when a small undetectable aggregate bonds with

Kinetics Analysis
While analysing the kinetics curves, three phases can be distinguished. The first one refers to nucleation, which happens on the molecular level, and does not cause the turbidity increase. The solution is clear, and mass growth is negligible. The second phase is visible as an intensive increase in absorbance value that can be correlated with silica skeleton formation via siloxane bonding. A plateau informs us of reaching the gel point where gel is in dynamic equilibrium-gel can be reorganised via coarsening and Ostwald ripening.
Within the model, the first phase (nucleation) cannot be observed in our numerical system-particles in CA are already formed, and their aggregation leads directly to the second phase. The third phase, regarded as a plateau, is when mostly aggregate clusters that were already detectable impact the structure formation process and morphology, but does not translate into an increase in the mass of the condensing product.
A summary of these theoretical considerations and assumptions is gathered in Table 2. We identified condensation mechanisms in each of the distinguished phases, according to Figure 5. During the first phase, the mass of the condensing product does not (or barely) increase. It means that the hydrolysed precursor molecules' condensation leads to nuclei formation, but they are too small to be detected by UV-Vis spectrophotometer. During the second phase, we have intensive mass growth. It means that large aggregates are becoming rapidly detectable. A mechanism that would explain such phenomena is bonding two undetectable aggregates into a bigger one which is immediately registered. Another option (however, less rapid) is when a small undetectable aggregate bonds with a detectable one. It is less rapid because the total detected mass would increase by one unit, but when two bigger (but not yet detectable) aggregates suddenly merge into a bigger one, the detected mass thus increases in step response. During the third phase, the mass is barely increasing because most of the clusters are big enough to be seen, and when they bond with each other, it does not impact mass growth, only structure formation and morphology.
The above qualitative considerations were quantified to compare and precisely analyse kinetics. The obtained results were interpolated by the spline method and differentiated. As described in our previous research [45], condensing skeleton mass is directly proportional to the absorbance value. Thus, we use a normalised absorbance value (A(t)/A max ) and identify its changes as condensation progress in relation to the final aerogel mass (A max ). The red line in Figure 6 is the first derivative (dm/dt), and it gives information about condensation rate evolution during the whole process. It also provides a specific and easy-to-determine point-the time when the maximum reaction rate is achieved, marked on dm/dt (blue line) as a grey dot. By making a tangent line to this point, it is possible designate the assumed time of the nucleation phase t 1 . Another point that we can quantify is the gel point. In our previous works, we observed that for particle-aggregate structures (nucleation and growth mechanism), the gel point t g was approximately when the normalised absorbance was equal to 99%. We assumed this point as the beginning of the third phase-the plateau. The period between these two points is the time of the second phase t 2 . a detectable one. It is less rapid because the total detected mass would increase by one unit, but when two bigger (but not yet detectable) aggregates suddenly merge into a bigger one, the detected mass thus increases in step response. During the third phase, the mass is barely increasing because most of the clusters are big enough to be seen, and when they bond with each other, it does not impact mass growth, only structure formation and morphology. The above qualitative considerations were quantified to compare and precisely analyse kinetics. The obtained results were interpolated by the spline method and differentiated. As described in our previous research [45], condensing skeleton mass is directly proportional to the absorbance value. Thus, we use a normalised absorbance value ( ⁄ ) and identify its changes as condensation progress in relation to the final aerogel mass ( ). The red line in Figure 6 is the first derivative ( / ), and it gives information about condensation rate evolution during the whole process. It also provides a specific and easy-to-determine point-the time when the maximum reaction rate is achieved, marked on / (blue line) as a grey dot. By making a tangent line to this point, it is possible designate the assumed time of the nucleation phase . Another point that we can quantify is the gel point. In our previous works, we observed that for particle-aggregate structures (nucleation and growth mechanism), the gel point was approximately when the normalised absorbance was equal to 99%. We assumed this point as the beginning of the third phase-the plateau. The period between these two points is the time of the second phase .

Condensation Kinetics Phases
Prevailing Mechanism of Structure Growth Experimental Numerical 1 nucleation mechanism A -(not considered in the model) 2 intensive mass growth (microscopic phase separation) phase mechanism B 3 cluster-cluster gelation (plateau) mechanism C The curves obtained experimentally (Figure 7) are sigmoidally shaped, with a short period of the first phase compared to the second one ( Figure 8). With the increase in the catalyst/precursor molar ratio, the gelation time decreases rapidly (A-D syntheses), then more slowly (E, F syntheses). When we look at the maximum value of the first derivative (condensation rate) (Figure 9), one can notice a drop for D and E syntheses. The dependence of dm/dt on time or reaction is presented in Figure 10. We can assume that three factors connected with dm/dt have an impact on the gelation time value: i.
The value of the maximum condensation rate ; ii. The time corresponding to the occurrence of the maximum rate; iii. The "span" of the first derivative's peak. The curves obtained experimentally (Figure 7) are sigmoidally shaped, with a short period of the first phase compared to the second one ( Figure 8). With the increase in the catalyst/precursor molar ratio, the gelation time decreases rapidly (A-D syntheses), then more slowly (E, F syntheses). When we look at the maximum value of the first derivative (condensation rate) (Figure 9), one can notice a drop for D and E syntheses. The dependence of dm/dt on time or reaction is presented in Figure 10. We can assume that three factors connected with dm/dt have an impact on the gelation time value: i.
The value of the maximum condensation rate dm dt max ; ii. The time corresponding to the occurrence of the maximum rate; iii. The "span" of the first derivative's peak.
Thus, even though the value of the first derivative is not constantly increasing, we observe that gelation time decreases with a higher contribution of catalyst in the system, presumably due to the imposition of the three (i-iii) factors. However, due to data interpolation accuracy, the first derivative seems to be the most reliable value that can be obtained numerically.             Kinetics curves obtained numerically for the mentioned 25 cases are presented in Appendix B. There is no first gelation phase, as nucleation is not considered during simulations. The dependence of the second phase duration on effective collision probability for variable concentrations of secondary particles in a system is presented in Figure 11. The dependencies have a power function character, shown as dotted lines. The values of the maximum condensation rate obtained for numerical simulations are presented in Figure 12. One can notice that the maximum condensation rate increases with increased particle concentration. It can be explained by referring to the mechanisms presented in Figure 6-the more particles we have, the more clusters will be created simultaneously, and a higher increase in reaction rate will be observed once the clusters reach the detectable threshold mass. However, there is no constant manner in the presented relation-the numerical kinetics curves are not as smooth and applicable for interpolation and differentiation as experimental curves obtained with low measuring time interval values. However, when looking globally at the numerical results (Appendix B, Figures 11 and 13), we see that reaction proceeds faster when: i.
The probability of an effective collision increases; ii. The concentration of secondary particles increases.
In Figure 13, numerical and experimental results of gelation time were plotted. Numerical t g values were designated as described in Section 2.3 as the point when the structure is 99% interconnected. Experimental t g values (triangles on the graph) were measured using the tilting test tube method but recalculated as a number of the time step so that it would be comparable with the numerical results (dots). To achieve that, the established dt values (Section 2.1) were used t g [dt] = t tilting g /dt.

Verification of Kinetics
Having in mind all of the adopted simplifications and assumptions within the model (discussed in the following Section 3), we attempted to compare the shape of the kinetics phase with a focus on the second phase of gelation, because it is the only phase that occurs for both the numerical and experimental system, and we can define the length of the phase based on the kinetic curve shape in a precise and repeatable way. We assumed that the duration of this phase is related to the maximum value of the condensation rate . The ratio / was calculated on the basis of experimental results and estimated for numerical data by extrapolating results for higher concentrations 6.2 % and interpolated for the probability of an effective collision per time step included in Table 3.
The results are presented in Figure 14. For the examined RLCA range, the ratio / depends linearly on the probability of an effective collision with convincing coefficients of determination ( = 0.87 for experimental and = 1 for numerical trendline). The formulas of the trendlines are presented in Table 3. The slope factors are very similar-the parameter for a numerical line is only 23% less than the experimental one. The difference for constant terms is higher-the value of numerical is positive, and the experimental is negative. Nevertheless, the disparities do not seem as significant when the linear character of the dependencies and the order of magnitude are the same.

Verification of Kinetics
Having in mind all of the adopted simplifications and assumptions within the model (discussed in the following Section 3), we attempted to compare the shape of the kinetics phase with a focus on the second phase of gelation, because it is the only phase that occurs for both the numerical and experimental system, and we can define the length of the phase based on the kinetic curve shape in a precise and repeatable way. We assumed that the duration of this phase is related to the maximum value of the condensation rate dm dt max . The ratio dm dt max /t 2 was calculated on the basis of experimental results and estimated for numerical data by extrapolating results for higher concentrations 6.2 % vol and interpolated for the probability of an effective collision per time step included in Table 3.
The results are presented in Figure 14. For the examined RLCA range, the ratio dm dt max /t 2 depends linearly on the probability of an effective collision with convincing coefficients of determination (R 2 = 0.87 for experimental and R 2 = 1 for numerical trendline). The formulas of the trendlines are presented in Table 3. The slope factors are very similar-the a parameter for a numerical line is only 23% less than the experimental one. The difference for constant terms is higher-the value of numerical b is positive, and the experimental is negative. Nevertheless, the disparities do not seem as significant when the linear character of the dependencies and the order of magnitude are the same. /t 2 ratio on probability of an effective collision for (a) experimental, (b) numerical approach.

Approach Established Dependence Coefficient of Determination
sumed that the duration of this phase is related to the maximum value of the condensation rate . The ratio / was calculated on the basis of experimental results and estimated for numerical data by extrapolating results for higher concentrations 6.2 % and interpolated for the probability of an effective collision per time step included in Table 3.
The results are presented in Figure 14. For the examined RLCA range, the ratio / depends linearly on the probability of an effective collision with convincing coefficients of determination ( = 0.87 for experimental and = 1 for numerical trendline). The formulas of the trendlines are presented in Table 3. The slope factors are very similar-the parameter for a numerical line is only 23% less than the experimental one. The difference for constant terms is higher-the value of numerical is positive, and the experimental is negative. Nevertheless, the disparities do not seem as significant when the linear character of the dependencies and the order of magnitude are the same. Figure 14. The dependence of (dm/dt)max/t2 ratio on probability of an effective collision for experimental and numerical approach. The dependence of (dm/dt) max /t 2 ratio on probability of an effective collision for experimental and numerical approach.

Structure
The porosity of alcogel samples (A-F syntheses) is presented in Figure 15 as a function of probability calculated based on the interpretation of the Arrhenius equation. The tendency can be approximated as a power function with a moderately fair R 2 value. The main observation is that porosity mainly increases along with probability value. One of the samples (E) could not be taken into account because of the exceptionally high measurement errors of calculation of probability (as presented in Figure 10). A similar observation regarding the porosity-probability relation was for numerical results (Figure 16), which can also be approximated by a power function. However, when the results are compared together (numerical vs. experimental), it can be noticed that alcogel porosity measured experimentally is much higher than that predicted via simulation results. The overall scheme presenting relations between particle concentration, probability of effective collisions and structure properties (morphology and porosity) is included in Figure 17. One can see how the structure morphology changes with model parameters (concentration and probability). Naturally, with a concentration increase, an aggregate becomes bigger. With increased probability, the structure tends to become more branched (DLCA structures) compared to more dense structures obtained for low probabilities (RLCA).

Structure
The porosity of alcogel samples (A-F syntheses) is presented in Figure 15 as a function of probability calculated based on the interpretation of the Arrhenius equation. The tendency can be approximated as a power function with a moderately fair value. The main observation is that porosity mainly increases along with probability value. One of the samples (E) could not be taken into account because of the exceptionally high measurement errors of calculation of probability (as presented in Figure 10). A similar observation regarding the porosity-probability relation was for numerical results (Figure 16), which can also be approximated by a power function. However, when the results are compared together (numerical vs. experimental), it can be noticed that alcogel porosity measured experimentally is much higher than that predicted via simulation results. The overall scheme presenting relations between particle concentration, probability of effective collisions and structure properties (morphology and porosity) is included in Figure  17. One can see how the structure morphology changes with model parameters (concentration and probability). Naturally, with a concentration increase, an aggregate becomes bigger. With increased probability, the structure tends to become more branched (DLCA structures) compared to more dense structures obtained for low probabilities (RLCA).
All of the obtained structures for the studied numerical cases are presented in Appendix A.

Discussion
The goal of the paper was to compare the experimental and numerical curves quantitatively. However, an awareness of several disparities in the experimental and numerical approaches should be raised. Within the numerical system, the initial size distribution of secondary particles is extremely monodisperse. The second phase of condensation proceeds, as it was subsequential to nucleation. Furthermore, the reaction is irreversible; no bond/interaction breaks within the already-created aggregate.
In the experimental system, the stages are not strictly subsequential in the same simple manner as in the cellular automaton. The borders can be indistinct for reactions that can occur parallelly or reversibly (hydrolysis, condensation) and for condensation phases (nucleation/primary particle formation, secondary particle aggregation). Massive aggregates can sediment. The experimental system is very complex, and the process cannot be observed in vivo from the very start (molecular level).
However, despite the assumed simplification at a certain level, results can be considered satisfactory and complementary to the current state of the art. The estimated values look consistent with the literature data for different organoalkoxysilane precursors (comparison is presented in Table 4).  Overall scheme presenting relations between particle concentration, probability of effective collisions and structure properties: morphology, fractal dimension and porosity.

Discussion
The goal of the paper was to compare the experimental and numerical curves quantitatively. However, an awareness of several disparities in the experimental and numerical approaches should be raised. Within the numerical system, the initial size distribution of secondary particles is extremely monodisperse. The second phase of condensation proceeds, as it was subsequential to nucleation. Furthermore, the reaction is irreversible; no bond/interaction breaks within the already-created aggregate.
In the experimental system, the stages are not strictly subsequential in the same simple manner as in the cellular automaton. The borders can be indistinct for reactions that can occur parallelly or reversibly (hydrolysis, condensation) and for condensation phases (nucleation/primary particle formation, secondary particle aggregation). Massive aggregates can sediment. The experimental system is very complex, and the process cannot be observed in vivo from the very start (molecular level).
However, despite the assumed simplification at a certain level, results can be considered satisfactory and complementary to the current state of the art. The estimated values look consistent with the literature data for different organoalkoxysilane precursors (comparison is presented in Table 4). All of the obtained structures for the studied numerical cases are presented in Appendix A.

Discussion
The goal of the paper was to compare the experimental and numerical curves quantitatively. However, an awareness of several disparities in the experimental and numerical approaches should be raised. Within the numerical system, the initial size distribution of secondary particles is extremely monodisperse. The second phase of condensation proceeds, as it was subsequential to nucleation. Furthermore, the reaction is irreversible; no bond/interaction breaks within the already-created aggregate.
In the experimental system, the stages are not strictly subsequential in the same simple manner as in the cellular automaton. The borders can be indistinct for reactions that can occur parallelly or reversibly (hydrolysis, condensation) and for condensation phases (nucleation/primary particle formation, secondary particle aggregation). Massive aggregates can sediment. The experimental system is very complex, and the process cannot be observed in vivo from the very start (molecular level).
However, despite the assumed simplification at a certain level, results can be considered satisfactory and complementary to the current state of the art. The estimated E values look consistent with the literature data for different organoalkoxysilane precursors (comparison is presented in Table 4). The research indicates that DLCA/RLCA model can be used for aggregation kinetics investigations. Furthermore, it was proved that experimentally measurable values for a certain synthesis, such as the time of gelation and alcogel porosity (or secondary particle diameter), can be translated into model parameters, such as the probability of an effective collision, particle concentration within a numerical system and time step of a simulation. According to our knowledge, such a "translation" of experimental values into numerical parameters was not yet described in the literature.
The similarity of the kinetics character shown in Figure 14 presents a very narrow range of probability values. According to our calculations, it is a very strong RLCA regime with a probability range 1.4 × 10 −4 , 3.1 × 10 −4 . The P range investigated within our cellular automaton studies = 1 × 10 −4 , 1 . Thus, discrepancies between tendencies observed for experimental and numerical systems can be observed. They were gathered and described in Table 5. Table 5. Similarities and discrepancies between tendencies observed for experimental and numerical systems. Decreases (most likely in an exponential manner; power and linear functions are also a good approximation) with an increase in both the n b /n p molar ratio and P.

Parameter
Decrease in power function manner with the increase in probability.

t g
Decreases with an increase in both n b /n p and P. Exponential and power trendlines are a good fit. However, we use an exponential approximation for the designation of E, A Decrease in power function of P.
Increases (but not constantly) with n b /n p or P. Seems to increase with concentration and P but in an unclear or constant manner.
The porosity comparison of structures obtained experimentally and numerically is not very satisfactory. The main tendency is the same (an increase in porosity along with the value of an effective collision probability); however, the values are underestimatedalcogel porosity measured experimentally is much higher than could be predicted based on simulation results. These aspects certainly should be improved in the future.

Experimental
The synthesis method was the two-step, acid-base sol-gel synthesis as described by Rao et al. [57]. The reaction mixture was prepared by mixing the proper volume of metyltrimethoxysilane (MTMS, purchased from Sigma-Aldrich, St. Louis, MO, USA) as a silica precursor, methanol (MeOH, Stanlab, Lublin, Poland) as a solvent and an aqueous solution of oxalic acid (concentration 0.01 M, Sigma-Aldrich) as an acidic catalyst to initiate the hydrolysis reaction. After a suitable time for this reaction to occur, an aqueous solution of ammonia (1 M, Eurochem BGD, Tarnów, Poland) was added to the mixture as the base catalyst of the gelation reaction. As a result, the alcogel was obtained. Six samples were chosen for the condensation kinetics investigation. All samples were prepared to have a constant MTMS/MeOH ratio (vol. 1:2) and a variable volume of aqueous catalyst solutions. Chemical compositions are presented in Table 6. The kinetics of gelation were investigated using measurements of UV-Vis light absorbance during the condensation process (scheme presented in Figure 18), utilising a Genesys 10 S UV-Vis Spectrophotometer Thermo Scientific with a wavelength λ = 633 [nm]. Changes in the absorbance of the mixture were found to be proportional to the mass of the formed alcogel, which was justified in our previous work [45].

Experimental
The synthesis method was the two-step, acid-base sol-gel synthesis as described by Rao et al. [57]. The reaction mixture was prepared by mixing the proper volume of metyltrimethoxysilane (MTMS, purchased from Sigma-Aldrich, St. Louis, MO, USA) as a silica precursor, methanol (MeOH, Stanlab, Lublin, Poland) as a solvent and an aqueous solution of oxalic acid (concentration 0.01 M, Sigma-Aldrich) as an acidic catalyst to initiate the hydrolysis reaction. After a suitable time for this reaction to occur, an aqueous solution of ammonia (1 M, Eurochem BGD, Tarnów, Poland) was added to the mixture as the base catalyst of the gelation reaction. As a result, the alcogel was obtained. Six samples were chosen for the condensation kinetics investigation. All samples were prepared to have a constant MTMS/MeOH ratio (vol. 1:2) and a variable volume of aqueous catalyst solutions. Chemical compositions are presented in Table 6.
The kinetics of gelation were investigated using measurements of UV-Vis light absorbance during the condensation process (scheme presented in Figure 18), utilising a Genesys 10 S UV-Vis Spectrophotometer Thermo Scientific with a wavelength = 633 [nm]. Changes in the absorbance of the mixture were found to be proportional to the mass of the formed alcogel, which was justified in our previous work [45].
Besides the changes in absorbance, the gelation time was measured to evaluate kinetic parameters in the Arrhenius equation. We used the basic, although easy, and commonly used tilting test-tube method [40,[58][59][60]. The measurement for each sample was performed three times and averaged for the three chosen temperatures of gelation (25, 35 and 45 °C). We used the analogical assumption, as in our previous work [61], that the gelation time is inversely proportional to the mean (or initial) reaction rate and, thus, to the reaction rate constant (Equations (1) and (2)).  Besides the changes in absorbance, the gelation time was measured to evaluate kinetic parameters in the Arrhenius equation. We used the basic, although easy, and commonly used tilting test-tube method [40,[58][59][60]. The measurement for each sample was performed three times and averaged for the three chosen temperatures of gelation (25,35 and 45 • C). We used the analogical assumption, as in our previous work [61], that the gelation time is inversely proportional to the mean (or initial) reaction rate and, thus, to the reaction rate constant (Equations (1) and (2)).

Numerical Model Description
The three-dimensional computing domain, divided into a regular grid of cells, represents investigated system. We investigated methyltrimethoxysilane (MTMS) gels with a "particle aggregates" type of structure, according to Nakanishi's division [14]. It is a hierarchical structure where precursor molecules undergo poly-condensation reactions and form primary particles that subsequently aggregate into secondary particles [62]. Thus, the model represents DLCA/RLCA mechanisms (cluster aggregation).
In our system, we assume that phase separation has already occurred due to nucleation and growth mechanism, and there is a certain amount of secondary particles. The initial conditions of a simulation (location of secondary particles) are pseudo-randomly generated according to a set concentration value c = N occupied /N cells ·100% (where N occupied is the number of secondary particles and N cells is number of cells within the computing domain). According to the cellular automaton definition, each cell is described by one of a finite number of states. There are only two possible states for each cell of the automaton-a cell can be occupied by a solid or a liquid phase.
Generally, in cellular automata systems with aggregation, we can distinguish three probability values determining the system dynamics: (i) the movement/diffusion probability P D , (ii) the reaction (effective collision) probability P, and (iii) the breaking probability P B [63].
The diffusion probability was determined based on the Stokes-Einstein equation for the diffusion of spherical particles ( D ∼ 1/r), which indicates that the diffusion coefficient is inversely proportional to the size of a particle. Thus, P D ∼ 1/r, where r is the radius of an aggregate (assumed to be a sphere), and N a is the number of secondary particles within the aggregate.
A collision occurs when two particles meet in adjacent cells (according to Moore's type of neighbourhood presented in Figure 19). However, in general, not every collision leads to the permanent attachment of the particle to an aggregate or two aggregates to each other. For this to take place, apart from the collision, a surface reaction followed by the formation of a "neck" must also occur. Such a collision is called as an effective one. The effective collision probability was determined based on the interpretation of the Arrhenius equation. To evaluate the Arrhenius constant and activation energy, it was assumed that the gelation time t g was inversely proportional to the reaction rate k, thus t g ∼ 1/k.  Thus, the time of gelation can be expressed in the form of Equation (2), where A is a constant dependent on the Arrhenius constant. This representation of the Arrhenius equation allows the establishment of A and E values [13,56], thus the reaction probability for a specific experimental synthesis. More information on the evaluation of model parameters is presented in the following sections. If the reaction between adjacent particles does indeed occur, these particles start to act as one aggregate (particles that belong to different aggregates may also react with one another).
The last parameter P B is the probability of breakage between adjacent, connected particles, although this aspect is not taken into a consideration within the presented model.

Kinetics Simulation
As mentioned in Section 1, DLCA/RLCA modelling was not originally intended to investigate the system's kinetics but to study the structure formation. In our case, it was the main objective to develop an algorithm that would be analogical and comparable with experimental measurements. Therefore, we focused on the mechanism of spectrophotometric measurement. During gelation, at first, the reaction mixture is transparent. As the reaction progresses, clusters of repeating units are formed, and when they reach a certain, not yet quantified size, they appear to be visible in the UV-Vis spectrum. The sample gradually becomes turbid, causing an increase in the sample's absorbance value.
In the presented numerical system, when an aggregate reached a detectable size (i.e., consisted of a set number of particles), its mass was added to the total sum of the condensed product and, thus, the condensation kinetics curve increased. An important aspect was defining a detectable size to correlate properly with experimental measurements. Three values were chosen to be tested: (i) a detectable size equal to 1% of a total number of secondary particles, (ii) equal to 2.5%, and (iii) equal to 4%. This approach for detectable size definition allows for the scaling of the parameter due to the concentration value or computing domain size. The kinetics curves registered for these three cases are presented in Figure 20A-C and an overall comparison is shown in Figure 20D. Referring to our previous studies [45], a gelation point for previously investigated synthesis was observed when absorbance reaches approximately 99% of its maximum value, and the gel point is usually located right after the inflexion of a curve (presented in Figure 6a and included in reference [45]). Considering all these facts, the detectable size was chosen as 1% of a system's total number of secondary particles. One should be aware of the arbitrary character of this parameter. However, it seems as a good initial assumption. The kinetics designation procedure continued until all particles were connected to form one big aggregate.

Structural Analysis
The obtained structures were studied in terms of their porosity , determined as the number of empty cells within a volume of a certain sphere (which is centred in the aggregate's centre of gravity) ( Figure 21). However, some inaccuracies had to be taken into

Structural Analysis
The obtained structures were studied in terms of their porosity ε, determined as the number of empty cells within a volume of a certain sphere (which is centred in the aggregate's centre of gravity) ( Figure 21). However, some inaccuracies had to be taken into account, including the fact that porosity is underestimated in the regions close to the centre of the aggregate (in particular, if we chose a distance smaller than the radius of a single particle, we would obtain a result of 0). At the same time, this value may be overestimated for greater distances from the centre, where the aggregate goes smoothly into empty space (there is no strict boundary, let alone a sphere-shaped boundary). Therefore, we needed to select a representative area in the aggregate under consideration. We chose the mean distance of each particle to the aggregate's centre as a reference point as it could be easily and precisely determined and because the local porosity value seemed approximately constant around this point. We assumed our representative area for porosity designation would be 80% r, 120% r . An exemplary plot presenting local porosity dependence on the distance from the aggregate's centre and the chosen representative area is presented in Figure 22.

Structural Analysis
The obtained structures were studied in terms of their porosity , determined as the number of empty cells within a volume of a certain sphere (which is centred in the aggregate's centre of gravity) ( Figure 21). However, some inaccuracies had to be taken into account, including the fact that porosity is underestimated in the regions close to the centre of the aggregate (in particular, if we chose a distance smaller than the radius of a single particle, we would obtain a result of 0). At the same time, this value may be overestimated for greater distances from the centre, where the aggregate goes smoothly into empty space (there is no strict boundary, let alone a sphere-shaped boundary). Therefore, we needed to select a representative area in the aggregate under consideration. We chose the mean distance of each particle to the aggregate's centre as a reference point as it could be easily and precisely determined and because the local porosity value seemed approximately constant around this point. We assumed our representative area for porosity designation would be 80% ̅ , 120% ̅ . An exemplary plot presenting local porosity dependence on the distance from the aggregate's centre and the chosen representative area is presented in Figure 22.

Conclusions
This paper successfully attempts to translate experimental synthesis conditions into CA model parameters in terms of kinetics curve comparisons.
Time of gelation measurements were used for and designation and subse- Figure 22. Scheme of porosity designation. Exemplary plot presenting local porosity dependence on distance from the aggregate's centre and the chosen representative area.

Conclusions
This paper successfully attempts to translate experimental synthesis conditions into CA model parameters in terms of kinetics curve comparisons.
Time of gelation measurements were used for A and E designation and subsequently for estimating the effective collision probability.
Alcogel porosity (or measurement of mean secondary particle diameter) is applicable for estimating particle concentration in the numerical system.
Gelation time can also be used to approximate the real value of a time step. According to our knowledge, such a "translation" of experimental values into numerical parameters had not yet been described in the literature, and we believe it can advance the current state of the art.
The comparison of kinetics curve shape is satisfactory. The ratio dm dt max /t 2 depends linearly on the probability of an effective collision for both experimental and numerical results. The slope factor difference is 10.1%, and the constant term difference is 22.8%. Nevertheless, the disparities seem acceptable, as the linear character of the dependences order of magnitude is the same.
Regarding the structure of aggregates, the main relation between alcogel porosity and probability is the same for both experimental and numerical results: an increase in porosity along with the effective collision probability value. However, we do not consider the quantitative comparison as satisfactory and it should be improved in the future.
The comparison of structures obtained experimentally and numerically is not as satisfactory. The main tendency is the same (an increase in porosity along with the effective collision probability value). However, the values are underestimated-alcogel porosity measured experimentally is much higher than could be predicted based on simulation results. These aspects certainly should be improved in the future.
Two possible improvements would minimise the simplification level of the model and should be discussed. The first one is the determination of gel point-it could be defined as a point when aggregates cannot rotate without hindrance. It would be analogical to the tilting test tube measurement method in some way, and it should be more precise than our simplified method. The second improvement would be considering the parallelism of primary and secondary particle formations. It could be completed using DLCA/RLCA aggregation and simultaneously adding new particles to the system.
To summarise, base-catalysed colloidal gel formation was successfully implemented as a CA model. The applicability of the numerical approach on the prediction of gels' condensation kinetics and their structure was evaluated, and it seems reliable within a strong RLCA regime. The developed model reflects the kinetics properly within the investigated chemical composition range. The structure analysis needs improvement. Ultimately, a relatively simple numerical model reflecting silica-based gel formation was developed and verified.