Roughness Effects of Crack Surfaces on the Elastic Moduli of Cracked Rocks

Crack surfaces are usually rough on various scales, and are sensitive to loading stresses and hence significantly affecting the mechanical properties of cracked rocks. We design a number of dry- and fluid-saturated numerical cracked samples to investigate the roughness influence of crack surfaces on the elastic stiffness. The fracture surface roughness is characterized by non-uniform fracture radii. We calculate the elastic moduli of cracked samples by finite-element simulation. Comparisons with the theoretical predictions by Gassmann and C&S (Ciz and Shapiro) (Ciz and Shapiro, Geophysics, 2007, 72(6), A75–A79) substitution equations demonstrate that the rough crack surfaces for both dry- and fluid-saturated samples can induce a stress concentration around the crack that reduces the elastic moduli and decreases the stiffness of rocks. For the fluid/solid-saturated cracks under the normal (shear) loading stresses, because the stress-concentration can induce shear (normal) strains around fracture, shear (bulk) modulus of the filling material will have contributions to the effective bulk (shear) modulus of rocks. The extra contribution, however, makes the Gassmann equation and C&S equation invalid.


INTRODUCTION
The characterization of fracture geometry is essential for a wide range of applications, such as geothermal production, hydrocarbon exploration, nuclear waste disposal, and CO 2 storage (Pruess, 2006;Hyman et al., 2015). The medium stiffness can be an effective tool for characterizing their surface geometry because of the sensitivity of medium stiffness to the fracture surface geometry, (e.g. Liu, 2005;Sevostianov and Kachanov, 2008;Gao and Gibson, 2012), which has been a strong interest of researchers.
Techniques aiming to detect fracture networks and surfaces in formations have attracted many attentions, (e.g. Kawahara, 1992;Sevostianov and Kachanov, 2008;Guo et al., 2018a). Fracture networks and surface geometry will influence the effective elastic moduli of rocks, which, in turn, can be used to detect the fracture networks and surface geometry features, (e.g. Zimmerman et al., 1986;Kachanov, 1993;Sevostianov and Kachanov, 2008). Researches on these issues can be traced back to around 1970s. Toksöz et al. (1976) take pores and flaws in porous rocks as oblate spheroids with varying aspect ratios, which tend to be closed under differential pressures, and the change of fracture density will influence the stiffness of rocks. Cheng and Toksöz. (1979) conduct a comprehensive investigation of stress-induced velocity variations using fractures with varying aspect ratios, and confirm that the deformation of fractures significantly influences elastic wave velocities. Then, more researchers have focused on the influence of fractures on the rock stiffness properties, (e.g. Zimmerman et al., 1986;Zimmerman, 1991;Shapiro, 2003;David and Zimmerman, 2012;Fu, 2017, Fu andFu, 2018). According to these studies, in conclusion, the open or closure of fractures will influence the fractured rock stiffness. The tectonic stress controls the closure or open of the fractures, and therefore changes the stiffness of the fractured rocks (Shapiro, 2003;Grechka and Kachanov, 2006;David and Zimmermann, 2012;Fu and Fu, 2018;Zong et al., 2020). However, many investigations are based on phenomenological models, relating the fracture geometry influence to its stiffness by empirical equations. The parameters of the empirical equations are obtained by fitting experimental measurements, which do not really explain the influencing mechanism of the fracture geometry on the rock effective elastic moduli.
To approach the elastic properties of fractured rocks more theoretically, the compliance of ellipsoidal fractures becomes a major issue, (e.g. Eshelby, 1957;Budiansky and O'connell, 1976;Kachanov, 1980;Sevostianov and Kachanov, 2002a). Hudson (1981) ;Hudson. (1988) takes fractures as circular inclusions in an elastic host medium and formulates the fracture influence on the elastic properties. Kawahara et al., (e.g. Kawahara, 1992;Kawahara and Yamashita, 1992;Kawahara, 2011) further regard fractures as first-order perturbations on the host medium and propose a model of frequency-dependent elastic moduli for aligned slit fractured sample based on the Foldy approximation (Foldy, 1945). This method has been extended to poroelastic media (Galvin and Gurevich, 2007;Galvin and Gurevich, 2009;Song et al., 2017a;Song et al., 2017b;Fu et al., 2020), focusing on the fracture lengths influence on the elastic moduli. Guo et al. (2018d) and Song et al. (2019); Song et al. (2020c) further study the fracture thicknesses influence on the frequency-dependent elastic moduli for fluid-saturated fractured rocks. Recently, Song et al. (2020a), Song et al. (2020b) have further researched the frequency-dependent elastic moduli for fluid-saturated fractured rocks with rectangular cracks and compressible fluid. However, these researches focus on the analytical solutions with the assumption of aligned fractures of the same scale.
Because fractures are distributed randomly with complex networks and rough surface, numerical simulations have been widely used to calculate effective elastic moduli, for example, finite-difference method (FDM) simulation for the variation of effective elastic moduli vs. fracture density (Saenger and Shapiro, 2002;Saenger et al., 2004), and the scattering effect of fractures (Vlastos et al., 2003;Vlastos et al., 2006;Vlastos et al., 2007). Recently, poroelastic finite-element method (FEM) simulation has been used to investigate the frequency-dependent elastic properties of porous rocks with intersecting fractures (Rubino et al., 2008;Rubino et al., 2016), and has indicated that the attenuation of P-and S-waves becomes dominant in the presence of fluid diffusion in the connected fractures, (e.g. Quintal et al., 2014;Rubino et al., 2015). Zhu and Shao. (2017) have done a complete review of the models describing the fractured rock effective moduli. However, all the researches studying the fractured rock effective elastic moduli above assume the fractures are penny-shaped or slit, with smooth surface.
The fracture surfaces in real rocks are rough, (e.g., Sevostianov and Kachanov, 2002b;Zong et al., 2017, Zong et al., 2020. The surface roughness significantly affects the fracture stiffness. Although the fracture surfaces are complex for natural rocks, their characteristics can be described by random functions, (e.g. Brown and Scholz, 1985;Gao and Gibson, 2012). Greenwood and Williamson. (1966) investigate the normal contact deformation problem of a rough surface and a flat surface based on the Hertzian contact theory, by assuming the Gaussian and exponential distribution of the asperity heights. Walsh and Grosenbaugh. (1979) apply the assumptions and methods of Greenwood and Williamson (1966) to study the effective compressibility of the two contacted rough surfaces and conclude that the effective stiffness is linearly proportional with the pressure. After that, more attentions are paid to the elastic mechanisms of rough fracture surfaces, (e.g. Johnson et al., 1972;Brown and Scholz, 1985;Priest and Taylor, 2000). A comprehensive review on these researches is made by Johnson (1985). The fracture compliance has been calculated for several types of irregularities, (e.g. Gao and Gibson, 2012;Sevostianov and Kachanov, 2012). However, these researches are limited to a single fracture with rough surface. For real rocks, the fractures are existent in rocks as fracture cluster, (e.g. Saenger and Shapiro, 2002;Saenger et al., 2004). However, there are no detailed researches about the fracture cluster with rough surface influence on the effective elastic properties of the fractured rocks.
In this paper, we study the dependency of fractured rock effective elastic moduli on the fracture surface roughness. We design seven fractured numerical samples. Each sample contains a series of fractures with the same perturbation of the square of the radius of the fracture. Different samples have fractures with different perturbation of the square of the radius. The square of fracture radius is controlled by the normal distribution function (Brown and Scholz, 1985;Adler and Thovert, 1999). To focus on the perturbation of the square of fracture radius influence on the fractured medium, we maintain the expectance of the square of the fracture radius and the aspect ratios of the fractures as constants, and vary the standard square of the fracture radius deviation to change the roughness of the fracture. We use the perturbation of square of fracture radius to represent the fracture surface roughness. The FEM simulation is used to compute the effective elastic moduli of each sample. We also assume the orientation and distribution of fractures are uniform to eliminate their effects. The calculated elastic moduli can be correlated to different levels of surface roughness. Comparisons with the Gassmann and C&S substitution equations predictions are conducted to evaluate the influence of the fracture surfaces roughness.

METHODOLOGY
The rough fracture surfaces can be generated by random functions, (e.g. Gangi, 1978;Sevostianov and Kachanov, 2008;Fu et al., 2020). On the other hand, the ellipses are widely used to characterize the fractures, (e.g. Eshelby, 1957;Cheng, 1978;Hudson, 1981;Kachanov, 1992;Guo et al., 2018a). In this work, fractures are considered as ellipses with the square of radii controlled by random functions, (e.g., Eshelby, 1957;Budiansky and O'connell, 1976;Sevostianov and Kachanov, 2012), as formed by Eq. 1, where x and y are the coordinates of the point in 2D surface. a 1, b 0.4. The aspect ratio (α bR/aR b/a) of the fracture is 0.4. In this research, we focus on the influence of the perturbation of square of factures radius (R 2 ) on the effective elastic moduli of the sample. We maintain the mean area of the fractures in the sample as a constant value, and vary R 2 using the normal distribution function with the expectation of R 2 as 9 mm 2 . Seven numerical samples, containing fractures with standard R 2 deviation (perturbation of R 2 ) ranging from 0 to 2.1 mm 2 , are generated. The side length L of each sample is fixed at 60 mm. The area of each fracture is a constant as πab < R 2 >, equaling 11.304 mm 2 , where < > represents the expectation value. Twenty-five fractures are inserted in each sample. The total fracture area for each sample is 282.600 mm 2 . We assume there are no other pores in the host medium, and the porosity of the sample is the ratio between the total fracture area and the sample area, determined to be 7.85%. We consider the orientations of fractures are distributed uniformly. The uniform distribution of the centers of fractures are also assumed, so the rock samples are macroscopically isotropic. The samples with different perturbation of R 2 are shown in Figures 1A,B.
For the fractures with rough surfaces, among the parameters describing fracture surface roughness, (e.g. Joint Roughness Coefficient (JRC), Z 2 , and fractal dimension), Z 2 is one of the most widely used statistical parameters in surface roughness analysis, and is shown to correlate well with JRC. The equation calculating Z 2 is given by (Barton et al., 1985) Z 2 1 n n 1 L where n is the number of surface profiles, L is the profile length, and x i is the coordinate of the surface profile in the point i. z i is the  Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 626903 height of fracture surface at point i. If the value of x i-1 -x i is a constant, as C. Eq. 2 can be simplified as In this research, we use the radius in point i, R i , replacing z i , and Eq. 3 can be rewritten as Therefore, it can be obtained that According to Zhao et al. (2019), surface with higher value of Z 2 is rougher. We have calculated the relation between Z 2 and perturbation of R 2 as shown in Figure 2. There is a positive linear relation between the perturbation of R 2 and Z 2 . Therefore, the perturbation of R 2 can represent the roughness of the fracture surface.
To compute the effective elastic moduli of the fractured samples, we load the static homogeneous stresses on the samples, and solve the static elastic equation at each node by FEM simulation (Saxena and Mavko, 2014), with the static elastic equations as and the stress σ ij is expressed as where e ij is the strain tensor, given by and e e kk , where u i is the displacement of each node of the sample, K and μ are the bulk and shear moduli of each node in the sample, respectively. The elastic moduli of fractures are different from those of the host medium. The elastic moduli of the host medium are the same for the all samples. Combining equations from Eqs. 6-9, we can calculate the stress σ ij and strain e ij in each node of the samples, by FEM simulation. Through the stress distribution, we can observe the stress concentration directly. All the FEM numerical simulations are conducted using commercial software (COMSOL 5.5), which can mesh the fractures automatically, and satisfy the stable condition. This kind of commercial software has been used widely (Li et al., 2009). To maintain the accuracy of the simulation, the smallest scale of the mesh is 0.05 mm, as shown in Figure 3. It should be noted that, although the 3-D simulation is more suitable for modeling real rocks, according to Rubino et al. (2008); Rubino et al. (2015); Rubino et al. (2016), and Saxena and Mavko (2014), the accuracy of 2-D numerical simulation is acceptable to extract the effects of fracture surface roughness on effective elastic moduli. In this research, we conduct numerical simulation on 2-D model.
To compute the effective bulk moduli of isotropic samples, we load a homogenous confining pressure P, as presented in Figure 4A. By solving Eqs. 6-9, the strain e ij and e at each node are obtained, and the effective bulk modulus K e of each sample is Similarly, we obtain the effective shear modulus μ e from the shear tests, by loading the homogenous pure shear stress τ 12 , as shown in Figure 4B. The effective shear modulus μ e is expressed as The fractures are important channels and spaces for water/ kerogen diffusion and enrichment. Especially for the rocks with low porosity and permeability, (e.g., Gale et al., 2014;Guo et al., 2018c), fractures are usually saturated with water or kerogen. It is important to study the fracture surface roughness influence on fluid/solid substitution process of fractured samples. The elastic moduli of dry sample are computed by FEM simulation, and then, we calculate the elastic moduli of the water-and kerogensaturated samples by FEM simulation, Gassmann equation and C&S (Ciz and Shapiro) equation (e.g., Gassmann, 1951;Ciz and Shapiro, 2007), respectively.
Gassmann equation (Gassmann, 1951) is expressed as where K sat w and μ sat w are the bulk and shear moduli of the watersaturated sample, respectively. K dry and μ dry are the bulk and shear moduli of the dry sample calculated by FEM, respectively, in this research. The pore space modulus M is given as where α 1-K dry /K B is the Biot-Willis coefficient (Biot and Willis, 1957). K B is the bulk modulus of the host medium. K f is the bulk modulus of water. The C&S equation (Ciz and Shapiro, 2007) is summarized as where K sat k and μ sat k are the bulk and shear moduli of the kerogen saturated medium. K k and μ k are the bulk and shear moduli of kerogen. μ B is the shear modulus of the host medium. Because Gassmann equation and C&S equation are valid to predict the effective elastic moduli of homogenous saturated rocks, comparison between the FEM simulation results and the Gassmann equation/C&S equation results will verify whether the two equations are valid, and whether the samples are homogenous. The difference between FEM simulation results and the Gassmann equation/C&S equation results indicates the changing of the heterogeneity of the sample, (e.g. Brown and Korringa, 1975;Saxena and Mavko, 2014). The differences should be increasing with the increase of the heterogeneity of the samples.
According to Saxena and Mavko. (2014), both Gassmann equation and C&S equations are set up based on reciprocal theory. These two kinds of equations assume samples are homogenous (Mavko and Mukerji, 2013). When confining pressure is loaded on samples, there are only elastic energy and normal strain induced by the normal stress of the sample, and only the bulk modulus of the inclusion has contributions to the effective bulk modulus. It is the same for the effective shear modulus. However, as verified by Saxena and Mavko. (2014), when the sample is heterogenous (the heterogeneity is induced by the distribution or geometry of the inclusion), if the normal stress is loaded on the sample, it will induce shear stress around the inclusions, and shear strain will be generated. Therefore, the shear modulus of the inclusions will have contribution to the effective bulk modulus. It is also the same for the effective shear modulus. All the explanations above mean the distribution or geometry of the inclusion maybe will make the Gassmann equation and C&S equation invalid (Mavko and Mukerji, 2013).

RESULTS
The effective elastic moduli K e and μ e are obtained by FEM simulation. The host medium bulk modulus K B 29.78 GPa, the host medium shear modulus μ B 22.30 GPa. For the dry numerical samples, the elastic moduli of filling material are zero. In the first stage, effective elastic moduli of dry samples are computed. Then, the water-and kerogen-saturated sample are considered. The bulk modulus of water K f is 2.25 GPa, and the shear modulus is zero. The bulk modulus of kerogen K k is 2.9 GPa, and shear modulus of kerogen μ k is 2.70 GPa (Mavko et al., 2009). Before doing the calculation, we do verify that all the numerical samples maintain the isotropic assumption, first. We conduct the numerical tests as shown in Appendix A. According to the results in Appendix A, the macroscopic isotropic assumption of each sample is maintained. In addition, we should also test the consistency between FEM simulation and Gassmann equation/C&S equation for homogenous and isotropic medium. The numerical test sample is presented in Figure 5. To maintain homogeneity and isotropy of the sample, all fractures are perfect circles, with radius as 1.895 mm to maintain the fracture area as 11.304 mm 2 . The total porosity of the test sample is the same as the numerical samples in the study as 7.85%.
We use FEM simulation to calculate the elastic moduli of dry sample, first. Then, we use FEM simulation, Gassmann equation and C&S equation to calculate the elastic moduli of water-and kerogen-saturated samples, respectively. The results are given by Table 1. The numerical error between the substitution equations and numerical results are less than 0.4% for the elastic moduli of kerogen-saturated sample, and the shear modulus of water-saturated sample. The numerical error between Gassmann equation and numerical result is 1.5% for the bulk modulus of water-saturated sample. The error in the bulk modulus might be induced by the interaction between fractures (Guo et al., 2018a). Both Gassmann and C&S equations are valid for homogenous and isotropic and homogenous medium. When we compare the difference between the FEM simulation results and substitution equation, if the differences between the numerical and substitution equation results are larger than the numerical errors, it should be induced by the fracture surface roughness.

Elastic Moduli of Dry Cracked Sample
In this section, we calculate the elastic moduli of dry numerical samples. Through FEM simulation, we obtain the variation of bulk and shear moduli of dry samples. As presented in Figures  6A,B, both bulk and shear moduli of the dry samples are decreasing sharply with the increase of the perturbation of R 2 . Because porosity is a constant, and there are no fluid or solid filling the fractures, the decrease of the bulk and shear moduli are resulting solely from the increase of the perturbation of R 2 . The sample with higher fracture surface roughness has lower stiffness. We further calculate the bulk and shear moduli of watersaturated sample by FEM and Gassmann equation based on the elastic moduli of dry numerical sample calculated in the last subsection, respectively. Figures 7A,B show the bulk and shear moduli variations vs. perturbation of R 2 of the watersaturated samples for both Gassmann equation predictions and FEM results. The comparisons between the results of FEM simulation (circle) and the Gassmann equation (square) show that the bulk moduli calculated by FEM simulation are less than those from Gassmann equation with a constant value. This kind of difference may be caused by the interaction between the fractures, which is as also observed by Zhao et al. (2016) Zhao et al. (2020) have also used the stress amplification and stress shielding to explain the influence of the interaction. Because the relative difference between the two kinds of bulk moduli are at about 1.5%, similar to the numerical error, we do not focus on the bulk modulus variation of water-saturated sample. We mainly focus on the shear modulus. The difference in shear modulus is obvious, and is larger than numerical error. The shear moduli obtained from FEM simulation are higher than the value of Gassmann equation ( Figure 7B). As shown in Figure 7D, the rough fracture surfaces will induce normal stress around the fractures, when shear stress is loaded. Because the saturated fluid has bulk modulus, and the normal stress around fractures will interact with the fluid, the bulk modulus of the filling fluid will increase the effective shear modulus. However, for Gassmann equation, as shown in Eq. 13, the effective shear modulus is only related to the shear modulus of dry samples. Therefore, the effective shear moduli of water-saturated samples obtained by Gassmann equation are less than the values of FEM simulation. As plotted in Figure 7B, with the increase of the perturbation of R 2 , the differences of shear moduli between Gassmann equation and FEM simulation are increasing. This indicates the contribution of the bulk modulus of the filling fluid to the shear modulus is increasing. In this study, the porosity of each sample is maintained as a constant, and the orientations of the fractures are also uniformly randomly distributed (the probability of different orientations are the same). Only the fracture surface roughness can generate this difference. The fracture with rougher surface will cause more normal stress around the fractures. The contribution of the bulk modulus of the filling fluid will also be increased. Gassmann equation is valid for homogenous medium (Saxena and Mavko, 2014). However, the difference between the numerical simulation and Gassmann equation results and the stress distribution in Figures 7A-D verify that the ellipsoidal fractures with rough surfaces will induce the local stress heterogeneity and strain heterogeneity of the sample. The local stress heterogeneity and strain heterogeneity will lead to the failure of the Gassmann equation (Mavko and Mukerji, 2013;Saxena and Mavko, 2014).

Elastic Moduli of Samples Saturated With Kerogen
We continue to analyze the two sets of kerogens-saturated results predicted by FEM simulation and C&S equation (Ciz and Shapiro, 2007) in this subsection, respectively. In Figures  8A,B, both bulk and shear moduli are decreasing with the increasing perturbation of R 2 . The results of FEM simulation are higher than the results of C&S equation. C&S equation is only valid for homogenous sample. According to Eq. 15, the effective bulk modulus is only influenced by the bulk moduli of the host medium and the filling solid. However, for the fractured sample in this research, as shown in Figure 8C, there are shear stresses around the fractures, when the normal stress is loaded. The shear stress will interact with the inclusions, and the shear modulus of the inclusion will have contributions to the effective bulk modulus of the sample. Therefore, the effective bulk modulus calculated by FEM simulation is higher than the result of C&S equation. To analyze the variation of shear modulus as plotted in Figure 8B, we also calculate the normal stress around the fractures. According to the normal stress distribution in Figure 8D, the bulk modulus of kerogen will also increase the effective shear modulus, and generate higher shear modulus, as calculated by the FEM simulation. Additionally, compared with the elastic moduli calculated by C&S equation, the variation of the elastic moduli calculated by FEM simulation are more stable. One reason for this phenomenon is that, although the perturbation of R 2 will decrease the stiffness of dry samples, when we load normal (or shear) stress around the fractures, the fractures with rougher surfaces will induce more shear (or normal) stress around the fractures, and the shear (or bulk) modulus of the filling solid will also have contribution to the effective elastic moduli. C&S equation does not calculate the contribution of the shear (bulk) modulus of the inclusions. Therefore, the effective elastic moduli calculated by the FEM simulation is higher than that obtained by C&S equation, and the FEM simulated elastic moduli are more stable. Frontiers in Earth Science | www.frontiersin.org April 2021 | Volume 9 | Article 626903 8 As verified by Saxena and Mavko. (2014), only when samples are homogeneous and isotropic, C&S equation is valid (Eqs. 15,16). According to the results in Figures 8C,D, because of the stress concentration around the fractures induced by the rough fracture surface (Kubair and Bhanu-Chandar, 2008), the strain distributions of saturated samples are heterogenous, and the C&S equation is invalid.

DISCUSSIONS
From the simulation results of dry numerical samples, we have observed linear decrease of elastic moduli of dry samples vs. the perturbation of R 2 of the fractures, as presented in Figures 6A,B. When samples are saturated with water and kerogen, because of the stress concentration caused by the rough fracture surfaces, the shear and bulk moduli of the inclusions will have contributions to the effective bulk and shear moduli, respectively. The fracture with rougher surface will generate higher stress concentration, and the shear and bulk moduli of the inclusions will have more contributions to the effective bulk and shear moduli, which will decrease the reducing speed of the elastic moduli.
In this study, we have also calculated the stress distribution of the dry samples with different fracture surface roughness, as shown in Figures 9A-D. Figures 9A,B show normal and shear stress distribution for the dry rock sample with the fracture radius standard deviation as 0 mm, for confining pressure and shear stress loaded tests respectively. Because there is no fluid or solid in the sample, the elastic moduli of fractures are zero, and the stress inside the fracture is low. Because the fracture surface is smooth, the stresses concentration around the fracture is low. The stress distribution heterogeneity is also low, and the elastic moduli is higher, as shown in Figures 6A,B. With the increase in the perturbation of R 2 , the stress concentration around fractures become larger, as shown in Figures 9C,D (standard R 2 deviation is 2.1 mm 2 ). As a result, the heterogeneity of the stress distribution of the sample is increasing, and the elastic moduli become lower ( Figures 6A,B). In this research, we have maintained the porosity unchanged, and the orientations of fractures are uniform. The increase of the stress concentration is induced by the increase of fracture surface roughness (perturbation of R 2 ). It also can be known that the stress concentration induces the decrease of elastic moduli of dry samples.
According to the analysis above, fracture surface roughness controls the stress concentration around the fracture surface and the heterogeneity of the stress. The higher perturbation of R 2 of fracture will induce higher stress concentration. The stress will concentrate around the peak of the fractures with large ratio between its length and its thickness (Kubair and Bhanu-Chandar, 2008). When the perturbation of R 2 (roughness) of the fracture is higher, it will generate more peaks with larger ratio, as well as higher stress concentration. Stress concentration will induce strain concentration around the fracture, and increase the strain in local region. When the loaded stress is the same, (e.g. Kubair and Bhanu-Chandar, 2008;Fu and Fu, 2018), the total strain will increase, because of the strain concentration. Therefore, the sample, with larger stress concentration, will be softer, and the elastic moduli will decrease.
According to the results in Figures 7, 8, the elastic moduli variation of saturated samples predicted by Gassmann and C&S equations are more drastic. According to Gassmann equation and C&S equation, the effective bulk moduli of samples saturated with fluid or solid are only related to the bulk moduli of the inclusions and host medium, and the effective shear moduli are only related to the shear moduli of the inclusions and host medium. However, because of the surface roughness of the ellipsoidal fractures, when stress is loaded on the sample, there is stress concentration around the fractures, the shear and bulk moduli of the filling medium will also have contributions to the effective bulk and shear moduli, respectively. The fractures with higher perturbation of R 2 will induce higher stress concentration, and the filling medium will have more contribution to the effective bulk and shear moduli. As a result, the elastic moduli obtained by FEM simulation is larger than that predicted by Gassmann equation and C&S equation, and the variation of the elastic moduli obtained by FEM are the more stable.
It should be noted that, the bulk moduli of the watersaturated samples obtained by FEM simulation is less than that obtained by Gassmann equation with a constant value. The difference may be caused by the interaction between fractures. This kind of difference are also observed in previous studies, (e.g. Guo et al., 2018a;Guo et al., 2018b;Cao et al., 2019;Cao et al., 2020).
According to the works of Zhao et al. (2016); Zhao et al. (2020), the interaction of the cracks (stress amplification) will induce the inconsistency between the substitution equations and numerical simulation in the elastic moduli of the sample. Comparing the stress distribution as shown in Figures 9A-D, the crack interaction is stronger for the cracks with rougher surface. The difference between the substitution equation results and the numerical simulation is also increasing with the crack surface roughness as shown in Figure 7B and Figures 8A,B. Therefore, there the difference between the substitution equation results and the numerical simulation is also induced by the interaction between the interaction between the cracks, and the rougher crack surface will induce higher interaction.
In this research, we study the substitution process in the limitation of linear elastic problem, and we do not consider the fracture surface slipping and the intersection between factures, when stress is loaded. Although fracture surface slipping and the intersection between factures will change fracture stiffness a lot, it is a kind of nonlinear elastic/ hypoplasticity problem (Jiang and Liu, 2008;Khidas and Jia, 2010;Jia et al., 2011). Generally speaking, both fracture surface slipping and the intersection between factures can increase the facture porosity and change the connection of the fractures (Shapiro 2003;Guo et al., 2018c), and they will change the structure and decrease the stiffness of the dry sample. In this research, we focus on the influence of fracture surface roughness on the substitution process, and we should maintain the fracture porosity and the stiffness of the dry sample as a constant during the stress loaded. Therefore, the fracture surface slipping, and the intersection influence are out of the region of this research. However, it should be noted that because the fracture surface slipping and the intersection between factures will change the stiffness of the dry sample and connection of the fracture during the stress loading process, it will also influence the substitution process. In addition, to focus on the fracture surface roughness influence, we assume all fractures are isolated to each other. However, for the for fractures in the natural world, it is possible that the fractures are intersecting to each other. When fractures are intersecting to each other, the interaction between the fractures will increase, and the stiffness of both dry and saturated samples will decrease. All the problem about the slip of fracture surface and fracture intersection will be researched in the future.

CONCLUSIONS
In this paper, we have studied the dependency of fractured sample elastic moduli on the fracture surface roughness. We design seven fractured numerical samples. Each sample contains twenty-five fractures with the same surface roughness. Different samples have different fracture surface roughness. The square of fracture radius is controlled by the normal distribution law. The FEM simulation is used to compute the effective elastic moduli of each sample. The resulting elastic moduli can be correlated to different levels of fracture surface roughness. Comparisons between the FEM simulation and theoretical predictions by the Gassmann and C&S substitution equations demonstrate that the fracture surface roughness could induce stress concentration, and thus reduce the elastic moduli of samples. For the sample with fractures orientated radomy uniformly, the factures with higher rough surface will induce more stress concentration around the fractures, and the interaction between fractures will increase, making the Gassmann and C&S equation invalid.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors. In this research, we make the orientations of the fractures are unform, assuming the samples are isotropic. To verify the correction of this assumption, we conduct the numerical test as shown in Figure A-1.
We load normal stress P in the two boundaries of the sample, for the dry rocks, and the other boundaries are made to be fixed, and we calculate the strain e xx and e yy for the test in Figures  A-1A, 1B respectively, and the elastic moduli in the two tests are given by C 11 P e xx .
(A-1) and C 22 P e yy , (A-2) where C 11 and C 22 are the elastic moduli of the samples in the two directions as shown in Figures A-1A, 1B. Then we calculate the relative error between C 11 and C 22 by Eq. A3 as e |C 11 − C 22 | C 11 . (A-3) The variation of the relative error vs. perturbation of R 2 is shown in Figure A-2. According to the results in Figure A-2, the relative errors are all less than 5%. Therefore, the samples can be regarded as isotropic.