Surface-Related Kinetic Models for Anaerobic Digestion of Microcrystalline Cellulose: The Role of Particle Size

In this work, for modelling the anaerobic digestion of microcrystalline cellulose, two surface-related models based on cylindrical and spherical particles were developed and compared with the first-order kinetics model. A unique dataset consisting of particles with different sizes, the same crystallinity and polymerisation degree was used to validate the models. Both newly developed models outperformed the first-order kinetics model. Analysis of the kinetic constant data revealed that particle size is a key factor determining the anaerobic digestion kinetics of crystalline cellulose. Hence, crystalline cellulose particle size should be considered in the development and optimization of lignocellulose pre-treatment methods. Further research is necessary for the assessment of impact of the crystalline cellulose particle size and surface properties on the microbial cellulose hydrolysis rate.


Introduction
Lignocellulose materials are abundant substrates for biogas and renewable energy production [1][2][3]. A better understanding of the key factors in the degradation process is required to increase the use of these materials for biogas production. Lignocellulose is a unique, recalcitrant structure of material plant cell walls. Its complex structure hinders microbial degradation. Further, cellulose is water insoluble and forms highly stable hydrogen-bonded crystalline fibers. These characteristics prevents its efficient utilisation in the anaerobic digestion process [4][5][6][7][8][9].
It is necessary to better understand the relationships between lignocellulose composition, degree of degradation, and degradation kinetics [10,11]. Moreover, there is still little information regarding the relationship between the kinetics of lignocellulose degradation and the structure of cellulose itself [12], particularly regarding the impact of crystallinity and enzyme adsorption sites [13]. The crystalline structure of lignocellulose hinders its degradation, and thus, the hydrolysis of these insoluble compounds is slow. However, the mechanisms behind the reduced hydrolysis rate are not fully understood; the low hydrolysis rate may be attributed to enzyme-associated factors such as enzyme inactivation, jamming due to enzyme overfilling on the cellulose surface, and inhibition of the enzymes. In addition, substrate-related factors, such as multiphase cellulose composition, changes in the degree of polymerization, crystallinity, reduction in reactivity, and substrate availability, as well as other physical properties may also contribute to the reduced hydrolysis rate [14].
In addition, the hydrolysis rate depends on many strongly related factors, such as the type and structure of the substrate (e.g., surface area and crystallinity), cellulose activity, and reaction conditions [15]. Some studies have shown that a decrease in crystallinity is invariably accompanied by a decrease in other substrate characteristics such as particle size and an increase in the available surface area that results from the pre-treatment of the material (mainly achieved by grinding) [16,17].
Hydrolysis of lignocellulose is considered to be the rate-limiting step in the anaerobic digestion process [11,18,19]. It has also been shown that the hydrolysis rate is controlled by enzyme kinetics if the amount of the enzyme exceeds the amount of the adsorption sites. In this case, hydrolysis can be described by first-order kinetics that have been widely used to describe the hydrolysis kinetics of cellulose and other solid particulate materials [20][21][22][23][24]. It is clear from the previous studies [21,22] that all accessible surface sites were used by enzymes; however, the first-order kinetics model does not describe a decrease in the amount of adsorption sites related to the shape and size of particulate matter. Furthermore, during anaerobic digestion, cellulose is hydrolyzed by cellulosomes that are multienzyme complexes produced by cellulose-degrading anaerobic microbes and mediate cell attachment to the of particulate [25][26][27]. After attachment to cellulose particle, cellulosomes grow rapidly and maintain contact with bacteria [26]; thus, enzyme-particle collisions can be necklaced as process rate limiting factor. Therefore, the available particulate surface is a crucially important factor. Some surface-related models of particulate hydrolysis have been developed [28][29][30], but this problem requires further investigation, particularly regarding cellulosome activity. To the best of our knowledge, surface-related models for anaerobic digestion have not been developed or improved upon in recent works. In this paper, we present a newly developed surface-related approach for the modelling of particulate matter degradation during anaerobic digestion using microcrystalline cellulose (MCC) as the model substrate. Pure cellulose was used in experiment to avoid interferences with lignin such as competitive and non-competitive inhibition [31]. The goal of this study was to isolate the role of cellulose particle size, so in our data set polymerization degree and crystallinity are also constant across samples, despite the fact that it is well known that mentioned above factors have significant impact on cellulose degradation kinetics [3]. It was assumed that the cellulosome activity on the entire accessible MCC surface was constant; however, particle surface, related to particle radius, decreased; thus, change in kinetics is strictly surface related. Based on this assumption, we developed two surface-related models for the hydrolysis of cylindrical and spherical particles, respectively, and compared these models to the first-order kinetics model. A unique dataset comprising of data for the particles with the same crystallinity and degree of polymerisation but with different sizes was used to validate the model. Such a dataset has not been used in previous studies. This constitutes an additional novelty in this paper. Due to a lack of reliable lag-phase-duration data pertaining to cellulosome formation as well as for the simplest possible modelling of the hydrolysis process activation, all of the tested models incorporated lag-phase modelling with the assumption of sudden hydrolysis process activation after the lag time.  [32,33]. The set of sieves had dimensions of the opening screens in the following sequence from the bottom to the top: 20, 32, 45, 56, 100, 150, 212, and 300 µm. A pan was located at the bottom. The sieve set was selected based on the initial test such that approximately half of the material passed through the middle mesh size. The mass of a single sample of cellulose subjected to separation was 50 g. The duration of sieving was 300 s and was controlled using a stopwatch. The end-point was determined based on preliminary experiments. For time durations between 300 and 600 s, the mass on the pan did not change by 0.1%; therefore, for economic reasons, it was decided to set the duration as 300 s. Seven samples of particle fractions between 20 and 300 µm were obtained. Mean geometric size of particles were calculated as geometric mean of top and bottom sieve size of given pair of sieves.

Anaerobic Digestion
The samples were processed by anaerobic digestion. The inoculum was obtained from a mesophilic, agricultural biogas plant post-digester. Inoculum contained 3.1% ± 0.05% total solids and 63.2% ± 0.05% volatile solids. Each fraction (10 g) was weighed on a WPS 600/C electronic scale (Radwag, Radom, Poland) with an accuracy of 0.01 g and placed into 2 dm 3 bottles. Following this, inoculum (1800 mL) was added. Bottles with inoculum only (1800 mL) were also prepared as blank. The inoculum-to-substrate ratio in the reactor was 3.5:1 based on organic dry matter in order to eliminate a possible negative effect of a high substrate dose on the process [34,35]. The bottles were flushed with nitrogen and placed in a water bath at 37 • C, and the reactor contents were manually mixed once a day. The amount of gas was measured by the brine displacement method [36]. Displaced brine was weighed on electronic scales with an accuracy of 0.01 g, as frequently as is necessary to fully recognize the course of gas formation (once a day or every second day at the end of process) [37]. Batch trials were simultaneously conducted in triplicate, and the obtained biogas volume was converted to the volume at standard conditions (1013 hPa, 273 K). The anaerobic digestion process was conducted for 25 days. Anaerobic digestion of the particle fraction in the 45-56 µm range was repeated because of trial failure on the first approach. Inoculum from the same biogas plant was used, but it was collected at a different time. Test results for particle fraction in the 45-56 µm range are given herein; however, as they were not obtained under the same conditions, they are not included in the overall analysis.

SEM Imaging
Images of randomly chosen microcrystalline cellulose samples were obtained by scanning electron microscopy (SEM, FEI Quanta 200 with microanalyser model EDS and digital image record; FEI Company, Hillsboro, OR, USA; merged with Philips Electron Optics) at the Analytical Centre of WULS. To obtain SEM images, a few cellulose particles were taken from each fraction and placed into a vacuum chamber.

True Density Measurement
Particle density was measured via gas stereopycnometry (Quantachrome Instruments, Boynton Beach, FL, USA). A mass sample m of MCC was placed in the measuring cell with a capacity of V C = 100 cm 3 . For this state, the pressure p 1 of the helium gas in the measuring cell was recorded. Then, the valve through which the gas was directed to the reference cell with the volume V A was opened, and its pressure p 2 was measured [38]. The true density of the MCC was determined as where ρ is the true density (g·cm −3 ), m is the mass of the MCC sample (g), V C is the volume of the measuring cell (cm 3 ), V A is the reference cell volume (cm 3 ), p 1 is the pressure in the measuring cell (MPa), and p 2 is the reference cell pressure (MPa).

Measurement of the Polymerization Degree
The polymerization degree was measured according to the PN-90/E-04421 and PN-92/P-50101/01 standards. Solution viscosity was assessed using an Ubbelohde viscometer (Schott Instruments GmbH, Mainz, Germany) with a capillary diameter of 0.84 mm. The cellulose solution outflow time was recorded every 10 min for 2 h for each sample weighing approximately 0.4 g. The capillary temperature was 25 ± 0.1 • C, and the error of the limiting viscosity number did not exceed 2%. The polymerization degree was calculated using the Immergut equation, DP α = K[η], where DP is the average degree of viscosity polymerization (-), [η] is the weight average intrinsic viscosity (-), and K and α are the empirical constants equal to 1.65 g·cm −3 and 0.9, respectively [39]. To minimize the deleterious effect of the copper (II) ethylenediamine (CUEN) solution, the polymerization degree for zero retention time was calculated by extrapolation.

Crystallinity Index Measurement
The crystallinity of the MCC samples (approximately 0.5 g) was evaluated by X-ray diffractometry (Tur M62 with a horizontal goniometer HZG-4, Carl Zeiss AG, Jena, Germany) with a Cu Kα radiation source (wavelength λ = 1.5418 Å) at a voltage of 30 kV and current of 25 mA. The scanning scope was 2θ = 5-30 • with a step size of 0.04 • , and the impulse counting time was 3 s. The crystalline structure of MCC gives rise to distinct peaks at 15 • , 16.4 • , and 22.5 • . The cellulose crystallinity was calculated by the XRD deconvolution method [40,41]. From the peak deconvolution method, the amorphous peak (2θ) was predicted to be located at approximately 21 • .

Specific Surface Measurement
Specific surface area (SSA) was measured using the Brunauer-Emmett-Teller (BET) method [42]. Samples were degassed in vacuum for 24 h in 105 • C prior to measurement in order to remove water and other contaminations. Measurements were performed using a QUADRASORB SI surface area analyzer (Quantachrome Instruments, Boynton Beach, FL, USA). Nitrogen gas was used as an adsorbate. All the necessary calculations were made using the 3 Flex analyzer software (v. 3.01, Micromeritics Instr., Norcross, GA, USA).

Water Retention Value Measurement
Water retention value (WRV) was measured according to the previously used protocol [43,44]. A precisely weighted sample (MA 50/1R electronic scale, Radwag, Radom, Poland with an accuracy of 0.1 mg) of cellulose (1 g) was placed in a centrifuge tube containing deionised water (10 mL), mixed, and left for 24 h. Then, the samples were centrifuged (3000× g for 20 min), and the supernatant was decanted. The remaining wet MCC samples were weighted and dried in an oven for 10 h in 105 • C. Then, the water retention value (WRV) was calculated as where WRV is the water retention value (g H 2 O·g −1 ), W wet is the weight of the wet centrifuged sample (g), and W dried is the sample weight after drying (g).

The First-Order Kinetics Model
In this study, three types of kinetic models were compared. The first was the first-order kinetics model with a lag phase [45] given by where B(t) is the biogas production at time t (mL·g −1 ), BP is the potential ultimate biogas production (mL·g −1 ), t is the time (h), k is the first-order kinetics constant, and λ is the duration of the lag phase (h). Sudden process activation after the lag phase was assumed. Thus, the lag during the lag time kinetics constant k is 0; the model is first deactivated, and then activated when t-λ is higher than 0. For the first-order kinetics model, the half hydrolysis time is expressed as where ln() denotes natural logarithm.
As mentioned above, this model focuses on the kinetics of substrate degradation with a kinetic constant k. It was found that the kinetic constant has different values for different substrate sizes even though this process occurs under the same pH and temperature conditions [21].
The SEM micrographs are presented in Figure 1. The SEM analysis of the shape of the MCC particles indicates that particle shape depends on the size of the mixture fraction separated by the sieves. Omitting surface irregularities and crystal roughness, it was assumed that the finer particles in the size range of 20-45 µm have a cylindrical shape and can be approximated as cylinders. In contrast, the particles with the sizes in the range of 56-212 µm form agglomerates in which the particle geometry can be approximated as spherical. These observations inspired the development of two mathematical models in which the access of enzymes to particles depends on the shape of the latter. Mathematical descriptions of the kinetics of the surface-related hydrolysis for the cylindrical and spherical particle shapes are presented below.
Materials 2021, 14, x FOR PEER REVIEW 5 of 18 As mentioned above, this model focuses on the kinetics of substrate degradation with a kinetic constant k. It was found that the kinetic constant has different values for different substrate sizes even though this process occurs under the same pH and temperature conditions [21].
The SEM micrographs are presented in Figure 1. The SEM analysis of the shape of the MCC particles indicates that particle shape depends on the size of the mixture fraction separated by the sieves. Omitting surface irregularities and crystal roughness, it was assumed that the finer particles in the size range of 20-45 μm have a cylindrical shape and can be approximated as cylinders. In contrast, the particles with the sizes in the range of 56-212 μm form agglomerates in which the particle geometry can be approximated as spherical. These observations inspired the development of two mathematical models in which the access of enzymes to particles depends on the shape of the latter. Mathematical descriptions of the kinetics of the surface-related hydrolysis for the cylindrical and spherical particle shapes are presented below.

Derivation of Surface-Related Models
A basic surface-related model equation was previously presented in the literature [46] as

Derivation of Surface-Related Models
A basic surface-related model equation was previously presented in the literature [46] as where M is the mass of the substrate (g), t is the time (h), S is the surface susceptible to hydrolysis (cm 2 ), and k s is the surface-related hydrolysis constant (g·cm −2 ·h −1 ); however, detail derivation of model assumptions was not presented there.
Step by step derivation is as follows.
It is possible to rewrite the above equation in terms of the released volume of the substrate, obtaining where V is the substrate volume (cm 3 ), and ρ is the substrate density (g·cm −3 ). Considering the degradation of spherical particles, the change in the volume can be expressed as where R is the particle radius (cm). We note that Thus, This simplifies to Thus, for spherical particles, it is possible to construct a model assuming an isotropic substrate and a constant decrease in the radius. It can be concluded that the rate of reaction is proportional to the particle surface area; thus, the decrease in the particle radius is constant. It was assumed that for cellulose, only the external layer of material is hydrolyzed at a given time. This approach is in line with the approach used for modelling the surfacerelated degradation of starch [46]. The present assumption for spherical particles is also valid for cylindrical particles with the assumption of radial particle degradation. Models sharing presented assumptions, adapted for biogas production curves, were not presented in literature before. In next sections full derivation of models adapted for biogas production curves are presented.

Model for Cylindrical Particles
Based on the presented assumptions for degradation kinetics, the model for cylindrical particles is expressed as where B(t) is the cumulative biogas yield at time t (mL·g −1 ), BP is the ultimate biogas production (mL·g −1 ), n is the number of particles, L is the average particle length (cm), and R is the initial average particle radius (cm). In this model, measured particle radius and length can be introduced; however, knowledge about particles number is then necessary.
The introduced lag phase model takes the following form The model shares lag-phase assumptions with the first-order kinetics model. If biogas or biomethane production curves are available without any additional information about the substrate specific surface or density, the only viable approach for modelling using surface-related models is relative modelling conducted with the assumption that the volume of the particle integrates to one. Additionally, assuming particles homogeneity, density is equal across particles. This allows the omission of the substrate density and the number of particles in the models. This makes models more clear in interpretation. Therefore, the model for the relative modelling of degradation for cylindrical particles is expressed as and where r s is the kinetic constant for relative modelling that considers the decrease in the particle radius. The unit of r s is cm·h −1 . It can be assumed that the total particle volume is equal to 1, ensuring that the most common structure of empirical equations is used for biogas or methane cumulative production curve modelling. Here, the constant related to the gas production potential is multiplied by an expression related to process dynamics. The dynamics-related expression integrates to 1 in its domain in nearly all cases [47]. It is assumed that the volume of a cylinder is equal to 1 Hence, Equation (14) can be transformed and simplified to where k is the first-order kinetic constant (h −1 ). Next, λ for lag-phase modelling is introduced, so that the final equation form is For Equation (17), the reciprocal of the constant k is the time required to complete the conversion after the lag phase. After this time, the constant k is assumed to be 0, as during the lag phase. The half decay time (t 0.5 ) is given by

Model for Spherical Particles
The model for spherical particles is based on the same assumptions as the previous model, with the exception that the particles are spherical; therefore, the basic equation is Introducing a lag phase, the model takes the following form For relative modelling, the model takes the following form and By analogizing to the model of cylindrical particles, it is possible to assume that volume of a sphere is equal to 1. Thus, Next, we include λ for lag-phase modelling; therefore, the final equation form is The reciprocal of the constant k is the time required to complete the conversion after the lag phase. The half decay time (t 0.5 ) is given by

Model Fitting
All of the calculations were performed using the R package (version 3.5.1). The functions implemented in the "GenSA", "minpac.lm" libraries were used to model biogas production. The models were fitted with the weighted least-squares method. The weights were calculated as a local slope of the biogas production curve [48] where w i is the weight factor (mL·g −1 ), B i (mL·g −1 ) is the biogas production curve value at a given measurement time i, and B i-1 (mL·g −1 ) is the previously measured value of the biogas production curve. Two algorithms were used during the approximation process. The Levenberg-Marquardt ("nlsLM" function from "minpac.lm" library) algorithm was preceded by a generalised simulated annealing algorithm ("GenSA" function from "GenSA" library) to initiate parameter estimation [49]. The default algorithm parameters were used for both algorithms. The global relative error was used to evaluate the fit [50] where n is the number of measurements, and z and p are the actual values and values obtained from the model, respectively. Error was estimated for the whole curves and alternatively for the last day of the trial only, to emphasize error on biogas yield. Because all the models used in this study had the same number of parameters, additional informational criteria (such as the Akaike criterion) were not determined.

Statistical Analysis
The Kruskal-Wallis test ("kruskal.test" form "stats" library) was used to compare the whole biogas production curves. The values of cumulative biogas production, cellulose moisture, and true density were compared using analysis of variance ("aov" function from "stats" library). Levene's test for the homogeneity of variance ("leveneTest" function from "car" library) and the Shapiro-Wilk normality test ("shapiro.test" function from "stats" library) were used to check the validity of the analysis of variance assumptions. Differences were considered statistically significant at p ≤ 0.05 for the statistical test.

Material Characterisation
The SEM images of MCC are shown in Figure 1, and it can be concluded that this material is characterized by irregular shape and polydispersity [51]. However, there are significant differences between the shapes of the small particles and large particles, with the geometric mean particle sizes of 25-38 and 75-178 µm, respectively. These particle shapes resemble cylinders and spheres, respectively. The smallest particles of MCC were mostly independent, while large particles formed spontaneous agglomerates due to the aggregation of single, fine crystals. Individual MCC particles from kenaf bast and wood pulp [52] and from cellulose filter papers [53] were determined to be rod-shaped. The particle shape of Flocel MCC with a mean volume diameter of 153 µm was determined to be generally fibrous with a few plate-like structures [54]. The physicochemical properties of the material are presented in Table 1. Significant differences were found for the true density of MCC particles (p = 0.0008), but the values varied in a very narrow range of 1.581-1.603 g·mL −3 . There was no evidence regarding differences in moisture (p > 0.05). The MCC moisture content was in the narrow range of 4.89-5.20% w.b., and is comparable to the values of 3.96-5.06% w.b. obtained for other types of MCC (Avicel, Flocel, fine powder, Ranq, fibres from sisal) [54]. The values of sample crystallinity can be compared when the moisture contents of the samples do not differ significantly [55]. XRD showed comparable results for crystallinity in the range of 56-59% (without the 45-56 µm particle fraction). In addition, the results for the polymerisation degree were comparable in a range of 328-336 (Table 1). Four types of MCC, namely, Avicel, Flocel, fine powder, and fibres from sisal, were obtained and had similar values of 60% for the crystallinity index (calculated by the Segal formula) [54]. The lack of correlation (r = -0.028) between the crystallinity and particle size differs from other findings in the literature [56]. However, it should be emphasized that the split of the mixture into fractions is not caused by the grinding of the material during which particle crushing occurs [56]. In this study, the finer particles merged into spontaneous agglomerates, and during this process, the fine particles did not undergo substantial physical changes. This indicates that the polymerization degree and material crystallinity do not impact the sample kinetics comparison. This conclusion is in line with previous studies that indicated that crystallinity is only one of several parameters that should be taken into account when assessing the enzymatic rate of cellulose degradation [57]. Further, relatively small changes in the crystallinity index should not be correlated with the changes in cellulose digestibility [41]. The results of the specific surface and WRV measurements are presented in Table 2. Statistically important differences in WRV were observed (p < 0.05). While the WRV values for the samples varied in a narrow range of 2.5-3.4 g H 2 O·g −1 , the differences were distributed across the whole range of particle diameter. Similarly, specific surface results varied only in the 0.73-0.87 m 2 ·g −1 range (Table 2) and were distributed across a range of diameters. In Section 3.3., the results presented in this section are analyzed in detail along with modelled data described in Section 3.2.

Results of Anaerobic Digestion and Model Comparison
Biogas production potential values are summarized in Table 3. Obtained biogas yield allowed to conclude that hydrolysis was completed, and theoretical biogas production potential was recovered [37]. Based on the Kruskal-Wallis test, significant differences (p = 0.0105) were found between the biogas production curves. Therefore, the curves were modelled separately. Biogas production curves are depicted in Figure 2, along with the model approximations for each MCC fractional sample. However, analysis of variance did not provide statistically significant evidence for the differences in cumulative biogas production (p = 0.0622). A comparison of the results of these two tests indicated that the differences between the biogas production curves can be attributed to different process kinetics.
Comparison of lag-phase duration, half decay time, and complete decay time is presented in Table 4. All models approximated the biogas production curves well. The obtained values for the lag-phase duration, half hydrolysis time, and complete hydrolysis time are presented. Lag-phase duration can be considered as an effect of whole population or single cell growth [58]. As previously mentioned, during the batch tests, the inoculumto-substrate ratio was 3.5:1 (based on organic dry matter); thus, the system was saturated with microorganisms. Therefore, the observed lag phase cannot be related to microbial population growth. The apparent lag should be considered as the time required for the synthesis of cellulosomes by individual cells and subsequent hydrolysis start-up. For the particles with the size in the 45-56 µm range, the lag-phase duration and half decay time were the longest, which was consistent for all of the tested models (Table 4).   This indicates that the inoculum used for anaerobic digestion of this sample, which was collected at a different time and had a lower activity, cannot be directly compared to the rest of the data. For this reason, the results for the particles with sizes in the 45-56 µm range were excluded from the further analysis of the lag phase, decay time, and overall kinetics. The experimental and model values of lag-phase duration were consistent for the rest of the particles and on average were 2.93, 2.56, and 2.71 d for first-order, cylindrical, and spherical particle models, respectively. This confirms the assumption of sudden hydrolysis process activation after the lag phase. This approach enables the calculation of the complete process time by adding the lag-phase duration and the time for the completion of decay/hydrolysis. Only surface-related models enable the calculation of the time to complete hydrolysis.
Comparison of model performance has shown that process kinetics are surface related, which is in line with previous research on cellulose hydrolysis by cellulosomes [13]. The spherical particle model estimated a longer hydrolysis time than the cylindrical particle model, and the average times were 15.64 and 11.40 d (excluding particles with the sizes in the 45-56 µm range), respectively. The smallest and largest differences were found for the particle fractions in the 20-32 and 212-300 µm size ranges (excluding particles with the sizes in the 45-56 µm range), respectively, and were 3.61 and 4.72 d, respectively (Table 4). The half decay times were similar for all tested models, with the values in the narrow range of 2.43-2.98 d; thus, differences between the surface-related models were observed at the end of the process.
The results shown in Figure 2 and the data presented in Table 3 show that the firstorder kinetics model overestimates real biogas production, with the exception of the curve for the 20-32 µm particle fraction, and the first-order kinetics estimates do not fit into the area marked by the curve and its standard deviation (Figure 2), in contrast to both surface-related models. This demonstrates the better prediction of anaerobic fermentation by these two models; however, it is not possible to clearly distinguish if the cylindrical or spherical model is better. The first-order kinetics model also exhibited the largest global relative error (6.03%) and the end-of-trial relative error (8.94%) for the 212-300 µm particle fraction. In comparison, the largest end-of-trial errors were observed for the 20-32 µm particle fraction, and were 2.37% and 2.10% for the cylindrical and spherical particle models, respectively (Table 5).
In most cases, the end-of-trial relative errors were slightly lower for the spherical particle model in than those obtained by the cylindrical particle model; the obtained values were 1.51% and 1.64% (excluding the fraction with the particle sizes in the 45-56 µm range), respectively, and can be attributed to the improved lag-phase modelling by the spherical particle model.

Analysis of Relations between Kinetic Constants, WRV and SSA
In this study, cellulose crystallinity and polymerization degree remained constant across samples; thus, k values were influenced by available adsorption sites related to substrate surface. For better analysis of k values, we proposed an empirical equation, which is the opposite Michaelis-Menten saturation equation where k is measured kinetic constant value (h −1 ), k max is the maximum kinetic constant value (h −1 ), K s is the half saturation constant (µm), and d is the mean particle diameter (µm).
The obtained global relative errors were 2.08% and 1.71% for cylinder and spherical shaped particles model, respectively. The results obtained using this approximation are presented in Figure 3. range), respectively, and can be attributed to the improved lag-phase modelling by the spherical particle model.

Analysis of Relations between Kinetic Constants, WRV and SSA
In this study, cellulose crystallinity and polymerization degree remained constant across samples; thus, k values were influenced by available adsorption sites related to substrate surface. For better analysis of k values, we proposed an empirical equation, which is the opposite Michaelis-Menten saturation equation where k is measured kinetic constant value (h -1 ), kmax is the maximum kinetic constant value (h -1 ), Ks is the half saturation constant (μm), and d is the mean particle diameter (μm).
The obtained global relative errors were 2.08% and 1.71% for cylinder and spherical shaped particles model, respectively. The results obtained using this approximation are presented in Figure 3. The proposed model interpretation refers only to process initialization and is as follows: when particle diameter is ≈0, nearly all adsorption sites are exposed. Then, with the increase in the particle diameter, the initial exposure of the enzyme adsorption sites de-  The proposed model interpretation refers only to process initialization and is as follows: when particle diameter is ≈0, nearly all adsorption sites are exposed. Then, with the increase in the particle diameter, the initial exposure of the enzyme adsorption sites decreases, so that the process requires more time based on the model assumptions (lower k constant). In this particular case, it can be stated that the difference in the particle size was the sole reason for better exposure of adsorption sites. During pre-treatment, both cellulose crystallinity and particle size usually change [16,17]. This is also true for ball milling of MCC [59]. A previous study on cellulose sonication has indicated the possibility of substantial hydrolysis enhancement due to considerable particle size reduction (from 38 µm to < 0.40 µm) with only slight changes in crystallinity (12% decrease) and degree of polymerization (no change), thus proving that recrystallisation is not inevitable [60]. Another study showed a considerable particle size decrease (from 300-500 to 20 µm in length and from 10-20 to 2 µm in diameter) after 120 h of enzymatic hydrolysis, while cellulose crystallinity decreased by only 12% [12]. Additionally, 400-500 nm channels were observed on the particles' surfaces after 48 h. Moreover, acid hydrolysis can lead to production of cellulose nanoparticles (size reduction from 45-53 µm to 53.5 ± 5.3 nm) with only 1.84% decrease in crystallinity, but a substantial decrease in the polymerization degree (from 466 to 215) [61]. Enzymatic hydrolysis can also lead to similar results, obtaining a size reduction from 45-53 µm to 36.5 ± 1.9 nm, polymerization degree decrease from 466 to 293, and 9.87% decrease in crystallinity [61]. Hence, it can be concluded that enzymatic hydrolysis causes only a 10-12% decrease in cellulose crystallinity, and a substantial decrease in the particle size. These observations are in line with our research, proving that change in the particle size only can enhance cellulose hydrolysis kinetics. This is confirmed by the recent studies on cellulosome activity that concluded that cellulosomes are particularly active for the smallest cellulose crystals (length ≤~70 nm) [62]. Such small particles are hardly attacked by free enzymes, possibly due to enzyme jamming [62,63]. However, for anaerobic digestion of lignocellulose, the influence of the particle size can be observed for much larger particles (0.15-1.70 mm) [64]. This can be attributed to the modification of the lignocellulosic matrix rather than the cellulose itself.
The WRV methods can be applied in order to measure the adsorption surface area. Typically, larger WRVs are associated with better surface exposure [65], larger inner pore volume, and lower crystallinity [43,44]. In MCC hydrolysis, hydrophobicity can be considered to be the key factor [44], and hence lower WRVs can be interpreted to be more beneficial if the crystallinity remains constant, which is true for this study. Figure 4 shows a comparison of WRV and k constant data, and it is observed that WRV increases with increasing particle size. Differences in uncertainty measures used on Figures 4 and 5 are the effect of different calculation software used. WRV changes can be attributed to the higher pore volume. However, WRV is correlated with the particle size (r = 0.822), and thus the WRV data cannot be interpreted without the knowledge of the particle size; however, it changes in a relatively narrow range (2.52-3.39 g H 2 O·g −1 ).
The BET specific surface area is frequently used to characterize lignocellulosic materials [43,66]. In this study, absolute SSA values do not differ substantially across samples, while kinetic constants show greater absolute differences, especially for the smallest particles. In this study, the surface area ranges from 0.73 to 0.87 m 2 ·g −1 (Table 2). Surprisingly, the specific surface area of the 20-32 µm particles (0.80 m 2 ·g −1 ) is lower than that of the 32-45 µm particles (0.87 m 2 ·g −1 ). This can be attributed to the particle structure, because the 20-32 µm particles are rather individual particles, while the 32-45 µm particles are agglomerates of smaller particles (Figure 1). While some studies have indicated that SSA is the most important indicator of hydrolysis effectiveness [44], this conclusion was based on ball-milled cellulose samples, where SSA was altered simultaneously with other properties such as a substantial crystallinity change. Additionally, during hydrolysis by free enzymes, not only the surface area, but also the surface properties, such as roughness, are important and can cause enzyme jamming [67]. Comparison of SSA with k constants is presented in Figure 5. In this study, the SSA does not accurately describe the changes in process kinetics.
In particular, the 20-32 µm particles show an opposite trend, because SSA decreases while k constant increases. Recent studies show that spatially organized enzymes in cellulosome can adapt their shape to cellulose nanocrystals. The individual cellulosome surface was calculated as approximately 1500 nm 2 [62]. SSA measured using the BET method is due to the areas of mesopores that are defined as the pores with the widths in the 2-50 nm range [68]. These pores are likely to be smaller than the proposed cellulosome area; hence, according to presented data, it can be hypothesized that cellulosomes do not necessarily utilize the mesopore surface. Additionally, the formation of agglomerates leads to decreased adsorption sites accessibility [69].
Materials 2021, 14, x FOR PEER REVIEW 15 of 18 cellulosome area; hence, according to presented data, it can be hypothesized that cellulosomes do not necessarily utilize the mesopore surface. Additionally, the formation of agglomerates leads to decreased adsorption sites accessibility [69].

Conclusions
All the tested models provide a good approximation of the hydrolysis process. However, the surface-related models outperform the first-order kinetics model that overestimates the true biogas production potential. The simple models presented in this work can be used as part of extended models in the evaluation of lignocellulose pretreatment meth-  cellulosome area; hence, according to presented data, it can be hypothesized that cellulosomes do not necessarily utilize the mesopore surface. Additionally, the formation of agglomerates leads to decreased adsorption sites accessibility [69].

Conclusions
All the tested models provide a good approximation of the hydrolysis process. However, the surface-related models outperform the first-order kinetics model that overestimates the true biogas production potential. The simple models presented in this work can be used as part of extended models in the evaluation of lignocellulose pretreatment meth-

Conclusions
All the tested models provide a good approximation of the hydrolysis process. However, the surface-related models outperform the first-order kinetics model that overestimates the true biogas production potential. The simple models presented in this work can be used as part of extended models in the evaluation of lignocellulose pretreatment methods. It is clear that the generation of amorphous cellulose during pretreatment is desirable because of its lower free energy demand for hydrolysis [44]. However, the unique data set presented in this paper and the newly developed surface-related modelling approach revealed that particle size is a key factor determining the kinetics of crystalline hydrolysis. Hence, crystalline cellulose particle size should be an important target property in the development and optimization of lignocellulose pre-treatment methods. Further research for the evaluation of the impact of the crystalline cellulose particle size and surface properties on the microbial cellulose hydrolysis rate is required.