## 2 Materials and Methods
### Generalized Singular Approximation Method (GSA)  
The general formula for effective stiffness tensor \(\mathbf{C}^{*}\) evaluation in the GSA method is [28; 29; 36]:  
\[\mathbf{C}^{*}=\Big{\langle}\mathbf{C}(\mathbf{r})\big{[}\mathbf{I}-\mathbf{g }\mathbf{C}^{\prime}(\mathbf{r})\big{]}^{-1}\Big{\rangle}\Big{\langle}\big{[} \mathbf{I}-\mathbf{g}\mathbf{C}^{\prime}(\mathbf{r})\big{]}^{-1}\Big{\rangle} ^{-1}, \tag{1}\]  
where angle brackets mean averaging over the representative volume; \(\mathbf{C}(\mathbf{r})\) is the fourth-rank tensor of the elastic constants related to the calcite matrix or fluid filling the voids; \(\mathbf{I}\) is the unit fourth-rank tensor; \(\mathbf{g}\) is the Green function tensor, defined by the equation:  
\[g_{ijkl}=-\frac{1}{16\pi}\int_{0}^{2\pi}\int_{0}^{\pi}\Big{(}n_{kj}\Lambda_{il }^{-1}+n_{ki}\Lambda_{jl}^{-1}+n_{lj}\Lambda_{ik}^{-1}+n_{li}\Lambda_{jk}^{-1 }\Big{)}\text{sin}\,\theta\text{d}\theta\,\text{d}\varphi, \tag{2}\]  
where \(\theta\in[0,\pi]\), \(\varphi\in[0,2\pi)\) and \(\Lambda_{ij}=C_{ijkl}^{c}n_{k}n_{l}\), \(i,j,k,l=1,2,3\) with \(n_{1}=\sin\theta\cos\varphi/a_{1}\), \(n_{2}=\sin\theta\sin\varphi/a_{2}\), \(n_{3}=\cos\theta/a_{3}\); \(a_{1},a_{2},a_{3}\) are semi-axes of the ellipsoidal inclusions. The tensor \(\mathbf{C}^{\prime}(\mathbf{r})\equiv\mathbf{C}(\mathbf{r})-\mathbf{C}^{\prime}\) is the fluctuation of the stiffness tensor. The tensor \(\mathbf{C}^{\prime}\) is the elastic tensor of the comparison body whose parameters could be arbitrary [28]. All indices in Equation (2) vary from 1 to 3. If we set \(\mathbf{C}^{c}\) equal to the stiffness tensor of the medium with the effective properties \(\mathbf{C}^{c}\)\(=\mathbf{C}^{*}\), then it gives us the formulas of the self-consistent method (GSA-SC) [28; 37]. This equation represents one of several realizations of self-consistent methods. The main idea of this approach is based on two hypotheses: (1) every inclusion acts like an isolated one in the medium with the effective properties of the composite and (2) the field acting on every inclusion is the external field applied to the medium. Thus, both hypotheses reduce the problem of interaction between many inclusions to a one-particle problem. In general, the self-consistent method gives good predictions for polycrystals and is a bit less accurate in the case of matrix composites. If we consider porous space connectivity in the \(f\)-model, we express the comparison body properties as a linear combination of the matrix and fluid properties according to Equation (3) [29]:  
\[\mathbf{C}^{c}=(1-f)\mathbf{C}^{\text{M}}+f\mathbf{C}^{\text{fl}}, \tag{3}\]  
where \(\mathbf{C}^{\text{M}}\) and \(\mathbf{C}^{\text{fl}}\) are, respectively, the stiffness tensors of the mineral matrix (the stiffest component) and fluid filling the voids (the softest component); \(f\) is an empirical dimensionless parameter reflecting connectivity of the voids. [PERSON] [38] shows that, in the general case, the parameter \(f\) is a tensor and the following expressions are valid:  
\[\mathbf{f}=\alpha\varepsilon^{\text{fl}}\Big{(}\varepsilon^{\text{C}}\Big{)}^{ -1},\mathbf{I}-\mathbf{f}=\alpha\varepsilon^{\text{M}}\Big{(}\varepsilon^{ \text{C}}\Big{)}^{-1}, \tag{4}\]  
where \(\mathbf{f}\) is a \"connectivity tensor\"; \(\alpha\) is the poroelastic Biot-Willis parameter [39; 40; 41]; \(\varepsilon^{\text{fl}}\), \(\varepsilon^{\text{M}}\) and \(\varepsilon^{\text{C}}\), are, respectively, the strain fields in the fluid inclusions, matrix, and com  
Figure 1: The illustration of proposed models of porous-cracked limestone: (**a**) double porosity model (Type 1 model), (**b**) \(f\)-model (Type 2 model). parison body. Thus, the physical meaning of the \(f\)-parameter specifies the deformation field distribution in the medium when the heterogeneous body reaches an equilibrium state. If \(f=1\) in Equation (4), we get the lower Hashin-Shtrikman bound (\(\mathbf{C}^{\mathrm{c}}=\mathbf{C}^{\mathrm{fl}}\)).