Azimuth-, angle-and frequency-dependent seismic velocities of cracked rocks due to squirt ﬂow

. Understanding the properties of cracked rocks is of great importance in scenarios involving CO 2 geological sequestration, nuclear waste disposal, geothermal energy and hydrocarbon exploration and production. Developing non-invasive detecting and monitoring methods for such geological formations is crucial. Many studies show that seismic waves exhibit strong dispersion and attenuation across a broad frequency range due to ﬂuid ﬂow at the pore scale known as squirt ﬂow. Nevertheless, how and to what extent squirt ﬂow affects seismic waves is still a matter of investigation. To fully understand its angle-5 and frequency-dependent behavior for speciﬁc geometries appropriate numerical simulations are needed. We perform a three-dimensional numerical study of the ﬂuid-solid deformation at the pore scale based on coupled Lamé-Navier and Navier-Stokes linear quasistatic equations. We show that seismic wave velocities exhibit strong azimuth-, angle-and frequency-dependent behavior due to squirt ﬂow between interconnected cracks. Furthermore, the overall anisotropy of a medium mainly increases due to squirt ﬂow but in some speciﬁc planes the anisotropy can locally decrease. We analyze the Thomsen-type anisotropic 10 parameters and adopt another scalar parameter which can be used to measure the anisotropy strength of a model with any elastic symmetry. This work signiﬁcantly clariﬁes the impact of squirt ﬂow on seismic wave anisotropy in three dimensions and can potentially be used to improve the geophysical monitoring and surveying of ﬂuid-ﬁlled cracked porous zones in the subsurface.

Abstract. Understanding the properties of cracked rocks is of great importance in scenarios involving CO 2 geological sequestration, nuclear waste disposal, geothermal energy, and hydrocarbon exploration and production. Developing noninvasive detecting and monitoring methods for such geological formations is crucial. Many studies show that seismic waves exhibit strong dispersion and attenuation across a broad frequency range due to fluid flow at the pore scale known as squirt flow. Nevertheless, how and to what extent squirt flow affects seismic waves is still a matter of investigation. To fully understand its angle-and frequency-dependent behavior for specific geometries, appropriate numerical simulations are needed. We perform a three-dimensional numerical study of the fluid-solid deformation at the pore scale based on coupled Lamé-Navier and Navier-Stokes linear quasistatic equations. We show that seismic wave velocities exhibit strong azimuth-, angle-and frequency-dependent behavior due to squirt flow between interconnected cracks. Furthermore, the overall anisotropy of a medium mainly increases due to squirt flow, but in some specific planes the anisotropy can locally decrease. We analyze the Thomsentype anisotropic parameters and adopt another scalar parameter which can be used to measure the anisotropy strength of a model with any elastic symmetry. This work significantly clarifies the impact of squirt flow on seismic wave anisotropy in three dimensions and can potentially be used to improve the geophysical monitoring and surveying of fluidfilled cracked porous zones in the subsurface.

Introduction
Wave propagation is controlled by the effective rock properties. Wave velocity and attenuation can be estimated from seismic data in scenarios such as seismic exploration, seismology, borehole measurements and tomography. Rock physics could then be used to estimate different rock properties, such as mineral composition, elastic moduli, the presence of a fluid, and pore space connectivity (and hence permeability) from seismic measurements. Thus, investigation of how cracks and fluids affect seismic properties has many practical applications. In activities including nuclear waste disposal, CO 2 geological sequestration, hydrocarbon exploration and production, geothermal energy production, and seismotectonics, a quantification of the fluid content, porosity and permeability of rocks are of great interest. All these activities can benefit from rock physics studies, and that is why cracked rocks have been under intensive studies during the last decades.
Published by Copernicus Publications on behalf of the European Geosciences Union. 856 Y. Alkhimenkov et al.: Azimuth-, angle-and frequency-dependent seismic velocities Fluid flow due to a passing wave may happen at different scales: at the wavelength scale, at the mesoscopic scale and at the pore scale (Müller et al., 2010). Biot's theory (Biot, 1962) describes the so-called global flow at the wavelength scale, but its overall effect on a passing wave at seismic frequencies is usually much smaller than that of fluid flow at the mesoscopic and pore scales (Pride et al., 2004). The mesoscopic scale is that much larger than the pore-scale but smaller than the wavelength. At this scale, studies are performed in the framework of Biot theory, assuming heterogeneous rock properties. One can define fractures as discontinuities at the mesoscopic scale and cracks as discontinuities at the pore scale. There are several analytical and numerical studies on the effect of wave-induced fluid flow between mesoscopic fractures and a porous rock background and between interconnected fractures using Biot's equations (Brajanovski et al., 2005;Rubino et al., 2013;Quintal et al., 2014;Masson and Pride, 2014;Grab et al., 2017;Hunziker et al., 2018;Caspari et al., 2019) as well as on the comparison between the numerical and analytical results (Guo et al., 2017(Guo et al., , 2018. Experimental studies of synthetic rock samples showed the impact of fluid-saturated fractures on seismic velocities (Amalokwu et al., 2016;Tillotson et al., 2012Tillotson et al., , 2014. The resulting frequency-dependent anisotropy was analyzed by Carcione et al. (2013), Rubino et al. (2017) and Barbosa et al. (2017). The last two also considered fracture-to-fracture flow, in addition to fracture-to-background flow.
At the pore scale, a passing wave induces fluid pressure gradients which occur between interconnected cracks, as well as between cracks and stiffer pores. Such pressure gradients force fluid to move between different cracks and pores until the pore pressure equilibrates throughout the connected pore space. This phenomenon, known as squirt flow (Mavko and Nur, 1975), causes strong energy dissipation due to the viscosity of the fluid and the associated viscous friction. Several experimental studies confirmed the importance of squirt flow at different frequency ranges (Mayr and Burkhardt, 2006;Best et al., 2007;Adelinet et al., 2010;Mikhaltsevitch et al., 2015;Pimienta et al., 2015;Subramaniyan et al., 2015;Chapman et al., 2019). There are several analytical solutions for squirt flow (O'Connell and Budiansky, 1977;Dvorkin et al., 1995;Chapman et al., 2002;Sarout, 2009, 2011;, which are based on simplified pore geometries and many physical assumptions.
Dispersion and attenuation caused by squirt flow can be simulated numerically by solving the coupled fluid-solid deformation at the pore scale using Lamé-Navier and Navier-Stokes equations with appropriate boundary conditions and then calculating effective frequency-dependent viscoelastic properties. During the last decades, many studies used numerical methods to solve mechanical problems (Andrä et al., 2013a, b;Saxena and Mavko, 2016). Recently, some numerical studies appeared in the geophysical literature aiming to solve the coupled fluid-solid deformation and hence studying dispersion and attenuation caused by squirt flow (Zhang et al., 2010;Zhang and Toksöz, 2012;Quintal et al., 2016Quintal et al., , 2019Das et al., 2019;Alkhimenkov et al., 2020). Das et al. (2019) numerically simulated a fully coupled fluidsolid interaction at the pore scale for digital rock samples. They modeled the pore fluids as Newtonian fluids using the Navier-Stokes equation with appropriate coupling between both the solid and liquid phases, accounting for inertial effects. Quintal et al. (2016Quintal et al. ( , 2019 simplified the Navier-Stokes equations by neglecting the inertial term and hence used the linearized quasistatic Navier-Stokes equation. We numerically simulate squirt flow in three dimensions and calculate frequency-dependent effective stiffness moduli using the finite-element method to solve the quasistatic Lamé-Navier equations coupled to the linearized quasistatic Navier-Stokes equations (Quintal et al., 2016Alkhimenkov et al., 2020). We apply an oscillatory deformation to certain boundaries of the numerical model, and, assuming that the wavelength is much larger than the size of individual cracks, we calculate the volume-average stress and strain fields and the resulting effective stiffness moduli. Then, we calculate the associated azimuth-, angle-and frequencydependent seismic velocities by solving the Christoffel equation. The main goal of this study is to analyze seismic anisotropy due to squirt flow in three dimensions since the previous numerical studies of seismic anisotropy were performed only in two dimensions and in the framework of Biot's theory Barbosa et al., 2017). This paper is organized as follows. First, we briefly describe the numerical methodology. Then, we describe the numerical model and show the numerical results -frequencydependent effective stiffness moduli. After, by solving the Christoffel equation, we evaluate the angle-, azimuth-and frequency dependent velocities of the model. Lastly, we quantify the anisotropy strength of the models analyzing the conventional Thomsen-type anisotropy parameters and also by adopting another scalar parameter.

Numerical methodology
We consider that at the pore scale, a rock is composed of a solid material (grains) and a fluid-saturated pore space (cracks). The numerical methodology is described by Quintal et al. (2016Quintal et al. ( , 2019 and Alkhimenkov et al. (2020), and here we briefly outline the main equations. The solid phase is described as a linear isotropic elastic material for which the conservation of momentum is (e.g., Lifshitz, 1959b andNemat-Nasser andHori, 2013) where "∇·" denotes the divergence operator acting on the stress tensor σ . The infinitesimal stress-strain relation for an where I 2 is the second-order identity tensor, tr is the trace operator, "⊗" defines the tensor product, the superscript "T" corresponds to the transpose operator, u is the displacement vector, and K and µ are respectively the bulk and shear moduli. The fluid phase is described by the quasistatic linearized compressible Navier-Stokes momentum equation (Landau and Lifshitz, 1959a): where v is the particle velocity, p is the fluid pressure and η is the shear viscosity. Equation (3) is valid for the laminar flow of a Newtonian fluid. In the finite-element numerical solver, Eqs.
(2)-(3) are combined in the space-frequency domain where ij represents the components of the strain tensor ij = 0.5 u i,j + u j,i , e is the trace of the strain tensor, λ and µ are the Lame parameters, u i is the displacement in the ith direction, δ ij is the Kronecker delta, i is the imaginary unit, and ω is the angular frequency. In the domain representing a solid material, Eq. (4) reduces to Eq. (2) by setting the shear viscosity η to zero. In the domain representing compressible viscous fluid, Eq. (3) is recovered by setting the shear modulus µ to zero. The solid and fluid displacements are described by the same variable, and thus they are naturally coupled at the boundaries between subdomains (Quintal et al., 2016. In the simulations, the energy dissipation is caused only by fluid pressure diffusion, since inertial terms are neglected. The COMSOL Multiphysics partial differential equation module is used for implementing Eqs. (1) and (4) (displacement-stress formulation) in a weak form. Our numerical results can be fully reproduced by using any openaccess software which includes mesh generation and finiteelement implementation with a corresponding solver for a linear system of equations. The whole spatial domain is discretized using an unstructured mesh with tetrahedral elements. A direct PARDISO solver (Schenk and Gärtner, 2004) is used for solving the linear system of equations. Direct relaxation tests are performed to compute all components of the stiffness matrix (in Voigt notation) c ij . The basic idea of the direct relaxation tests is that a displacement boundary condition of the form u = 10 −6 × exp(iωt) is applied to a certain external wall of the model and in a certain direction, while at other walls of the model the displacements are set to zero or let free to change. In the direct tests that we perform, only one component of the stiffness matrix c ij can be directly calculated after one numerical simulation. A detailed description of the boundary conditions is given in Alkhimenkov et al. (2020). The initial conditions for displacements are set to zero. The resulting stress and strains are averaged over the spatial domain for each frequency. Then, the complex-valued c ii (ω) components (diagonal) are calculated for each frequency (in Voigt notation, no index summation): where · represents the volume averaging over the sample volume.
Note that usually the inverse quality factor is used as a measure of attenuation (O'connell and Budiansky, 1978).
In this study, we show the inverse quality factor for each component of the stiffness tensor, even though the ratio Im c ij (ω) /Re c ij (ω) does not represent attenuation of any corresponding wave mode for some components.

Numerical model
Two 3D numerical models are constructed, which consist of a pore space embedded into an elastic solid grain material (Fig. 1). The solid grain material is represented by a cuboid whose size is (0.24×0.24×0.24) m 3 . The pore space consists of two perpendicular cracks represented by thin cylinders of 0.002 m thickness, 0.1 m radius (i.e., the aspect ratio is thickness divided by diameter -0.01) and fully saturated with a liquid. In the first model, the two cracks are disconnected, while, in the second model, the two cracks are connected (cross sections in Fig. 1). The employed liquid properties are those of glycerol, and the grain material has properties of quartz (Table 1). A fine, regular mesh is used inside the crack to accurately account for dissipation, while in the grain material the mesh is coarser (Fig. 2). The total number of elements is 3.3×10 6 . The simulation is performed for 12 different frequencies  One crack embedded into an isotropic background induces a transverse isotropy (five independent components of the stiffness tensor, e.g., Mavko et al., 2009). If the crack is parallel to the x-y plane, then the symmetry is vertical and the medium exhibits vertical transverse isotropy -VTI symmetry. If the crack is parallel to the x-z plane, then the symmetry is horizontal and the medium exhibits horizontal transverse isotropy -HTI symmetry. If two cracks, perpendicular to each other, are embedded into an isotropic material and the crack compliances are different, then the medium exhibits orthorhombic symmetry (nine independent components of the stiffness tensor). If the crack compliances are the same, then the medium symmetry is tetragonal (six independent components of the stiffness tensor); some authors attribute this geometry to a special case of orthorhombic symmetry (e.g., Bakulin et al., 2000b), while tetragonal and orthorhombic symmetry classes are different. On the other hand, one can argue that an orthorhombic medium (created by two perpendicular sets of cracks) degenerates into a tetragonal medium if the crack compliances are the same.
The symmetry of the saturated numerical model with connected cracks is tetragonal ( Fig. 1), because the crack compliances are the same. Thus, there are only six independent components of the stiffness tensor. We will see that the symmetry of the saturated numerical model with disconnected cracks is orthorhombic, because one crack is stiffer than the other one due to its separation into two parts. However, the difference between c 22 and c 33 stiffness components is less then 0.3 %, thus the divergence from the tetragonal symmetry is negligible, and therefore this model is considered tetragonal as well.

Dry stiffness moduli
Let us first consider the geometry shown in Fig. 1 with a pore space filled with air (i.e., dry). We perform nine relaxation tests to calculate the full stiffness tensor for each of the two models with connected and disconnected cracks. The result-ing effective stiffness moduli for the model with connected cracks are (in Voigt notation) For the model with disconnected cracks, the effective stiffness moduli are (in Voigt notation) The effective stiffness moduli of the two models are different. Zero values are written if the value is below 0.0002 GPa (i.e., up to numerical precision). The c Con ij stiffness matrix precisely belongs to the tetragonal symmetry class, while the c Dis ij stiffness matrix has all diagonal components different from each other; thus, it represents the orthorhombic symmetry class. The largest difference between c Con ij and c Dis ij is in the c 22 component; i.e., c 22 = c Dis 22 − c Con 22 = 5.3 GPa. That is a significant difference and it is only due to the vertical crack separation. There are two different features which must be clearly separated. (1) The effect of crack intersection without changing the crack geometry on the effective elastic properties. In this case, the crack intersection is achieved by changing the spatial position of the cracks. Grechka and Kachanov (2006) studied numerically the effect of crack intersection without changing the crack geometry. They concluded that crack intersection has a very little impact on the effective elastic moduli. (2) The effect of the crack partition into two "halves" on the effective elastic properties. In this case, the partitioned crack has a long, thin contact area across the whole diameter ( Fig. 1). It is well known that the contact areas inside a crack significantly reduce crack compliance (Trofimov et al., 2017;Kachanov and Sevostianov, 2018;Markov et al., 2019;Lissa et al., 2019). Comparing Eqs. (7) and (8), we also observe that the thin contact area significantly reduces the crack compliance: the effective dry moduli of the model with disconnected cracks are much higher compared to the model with connected cracks. An intuitive explanation is the following: if crack surfaces have not been changed by changing the spatial position of the cracks in the volume, the effect of crack intersection is negligible (Grechka and Kachanov, 2006); if the crack surfaces have been changed -as we did in the present study by partitioning the vertical crack into two pieces (and introducing a thin additional contact area) -the effective elastic moduli would become much stiffer compared to the model where the crack surfaces have not been changed.

Fluid pressure fields
Here and later on we deal only with a liquid-saturated pore space. The liquid has properties of glycerol (Table 1). A direct P-wave modulus test is performed to calculate disper-sion and attenuation for the c 33 component. (A harmonic displacement is applied to the top wall of the model in the z direction, while the normal component of the displacement is set to zero on all the other walls.) Figure 3 shows snapshots of the fluid pressure P f in the cracks at three different frequen-860 Y. Alkhimenkov et al.: Azimuth-, angle-and frequency-dependent seismic velocities cies, in the vertical middle slice of the model (the y-z plane, red frame in Fig. 1a). For the model with connected cracks, at low frequencies, there is enough time for pressure equilibration between the cracks; thus, the pore pressure is uniform throughout the pore space (Fig. 3, LF (connected)). This is called the relaxed state. At intermediate frequencies, there is a large pressure gradient in the cracks, which corresponds to the maximum attenuation due to squirt flow between cracks (Fig. 3, Fc (connected)). At high frequencies, there is no time for fluid to move; hence, there is no fluid pressure equilibration between the vertical and horizontal cracks (Fig. 3, HF (connected)). This is called the unrelaxed state. Therefore, at high frequencies, the connected cracks behave as hydraulically isolated and the fluid highly stiffens the crack.
In the model with disconnected cracks, the fluid pressure in the cracks is the same in all three regimes, which corresponds to the unrelaxed state in the model with connected cracks. The unrelaxed state can be interpreted as the elastic limit because there is no fluid flow between the cracks, and the effective properties of the two models (connected and disconnected cracks) are the same, as will be shown in the next subsection. Figure 4 shows the numerical results for the complexvalued frequency-dependent components of the stiffness matrix c ij (ω) (in Voigt notation) for the models with connected and disconnected cracks filled with glycerol. In the model with connected cracks, the real part of the c ij component and the corresponding inverse quality factor (Eq. 6) curves show strong frequency-dependent behavior of the c 22 , c 33 and c 23 components ( Fig. 4a and c). The inverse quality factor and dispersion of the c 22 and c 33 components coincide because the geometrical properties of the two cracks are the same (Fig. 4a) and the model is symmetric. The c 11 component is nondispersive and exhibits zero attenuation. The dispersion of the c 44 , c 55 and c 66 components is negligible and these components also exhibit negligible attenuation (Fig. 4b). The c 12 and c 13 components are nondispersive, the c 23 component exhibits strong negative dispersion and a negative inverse quality factor peak is shifted towards high frequencies compared to that of the c 22 and c 33 components. A similar phenomenon has been reported by Guo et al. (2017) in the context of two-dimensional simulations. The c 23 component does not correspond to a wave mode alone; it is always used together with c 22 and/or c 33 components. Therefore, no wave will gain energy. This negative inverse quality factor sign for the c 23 component was also verified using Kramers-Kronig relations. In other words, different components of the stiffness tensor might have positive or negative values of the ratio Im(c 23 )/Re(c 23 ) but, when we calculate the velocity and the inverse quality factor of a wave, the cumulative effect of all c ij components must be physical and no negative attenuation will be observed.

Elastic moduli
Note that the width of the inverse quality factor peak (at half amplitude) for the components c 22 and c 33 has a 1.5 order of magnitude ( Fig. 4a and 4c). It means that attenuation and dispersion due to squirt flow play a significant role over a broad frequency range even for cracks with a single aspect ratio.
In the model with disconnected cracks, all components of the stiffness tensor c ij (ω) (Fig. 4a and 4c) are constant across the whole frequency range and exhibit zero inverse quality factor. Furthermore, all components are approximately equal to the high-frequency values of the model with connected cracks. This is expected in the unrelaxed state because the connected cracks behave as hydraulically isolated with respect to fluid flow. A very small discrepancy between the two models at high frequencies is associated with the vertical crack partition (two thin regions of pore space replaced with stiffer grain material). Figure 5 shows the P-wave (primary wave) phase velocity as a function of the phase angle of the numerical model with connected and disconnected cracks (Fig. 1), where the zero phase angle corresponds to the vertical wave propagation (along z axis). The P-and S-wave phase velocities are calculated by solving the Christoffel equation, which represents an eigenvalue problem relating the stiffness components c ij , the phase velocities of plane waves that propagate in the medium and the polarization of the waves (Fedorov, 1968;Tsvankin, 2012). Considering the plane Y -Z, the P-wave velocity is the same for phase angles of 0 and 90 • ; it changes with frequency only for phase angles between 0 and 90 • and is maximal in the high-frequency limit at phase angle of θ = 90(±90) • (Fig. 5a). Furthermore, in the high-frequency limit the Pwave phase velocity coincides for the models with connected and disconnected cracks. As frequency decreases, the P-wave velocity decreases, and at 10 4 Hz the P-wave velocity is almost angle independent (yellow curve, Fig. 5a). It is interesting that this "local" isotropy corresponds to the maximum attenuation of the c 22 and c 33 components (Fig. 4). As frequency further decreases, the P-wave velocity decreases and stays nearly unchanged for the frequencies below 10 3.5 Hz.

P-and S-wave velocities
In the X-Z plane, the P-wave phase velocity is the same for the models with connected and disconnected cracks in the high-frequency limit (Fig. 5b). For the model with connected cracks, as frequency decreases, the P-wave velocity decreases, reaching its minimum at low frequencies (10 1 -10 3.5 Hz). Figure 3. Snapshots of the fluid pressure P f in the cracks at three different frequencies: LF -the low-frequency limit (corresponds to 10 1 Hz, relaxed state), Fc -intermediate frequency snapshot (corresponds to 10 4 Hz, close to the characteristic frequency) and HF -the high-frequency limit (corresponds to 10 6.5 Hz, unrelaxed state). Figures 6-7 show the quasi-shear (SV) and the pure shear (SH) phase velocities as functions of the phase angle of the numerical models with connected and disconnected cracks (Fig. 1). The SV wave velocity is strongly frequency dependent in both the X-Z and Y -Z planes. The SH wave exhibits some frequency-dependent behavior in the X-Z plane and is angle and frequency independent in the Y -Z plane. It is interesting that the SV waves in two different planes have different velocities at 0 and 90 • phase angles, which is due to their different wave polarization. The SV wave in the Y -Z plane has the same polarization as the SH wave in the X-Z plane; their velocities are equivalent at the 0 and 90 • phase angles. The same conclusion is valid for the SV wave in the X-Z plane and the SH wave in the Y -Z plane. A slight discrepancy (around 0.5 %) between the SV wave velocities for the disconnected crack model (Fig. 6, dashed red line) and the high-frequency velocity for the connected crack model (Fig. 6, green line) at phase angles of 0, 90 and 180 • is due to the crack separation.
Due to the symmetry of the model, the behaviors of the P-, SV-, SH-wave phase velocities in the X-Z and X-Y planes are identical; thus, the results in the X-Y plane are not shown here.

Quantitative analysis of the frequency-dependent anisotropy
First, we quantify the Thomsen-type anisotropic parameters (Thomsen, 1986) for orthorhombic media (Tsvankin, 1997;Bakulin et al., 2000a). Then, we quantify the universal elas-tic anisotropy index (Ranganathan and Ostoja-Starzewski, 2008) and the two parameters which define the anisotropy strength in bulk and shear modes. All these anisotropy measures highlight different frequency-dependent features of the models. Our results shown in Fig. 4 (frequency-dependent elastic moduli) are used as input to quantify these anisotropy measures.

Thomsen-type parameters
Thomsen-type anisotropic parameters ( , δ, γ ) describe the P-wave anisotropy -, the shape of the P-wave phase velocity at different phase angles -δ and the S-wave anisotropy γ : each set of three parameters corresponds to one plane. Thus, for our model symmetry, there are two different planes -Y -Z and X-Z (because the X-Z plane is equivalent to the X-Y plane). In this study, we refer to Thomsen parameters | |, |δ|, |γ | ∈ [0, 0.1] as to weak elastic anisotropy (| · | corresponds to the absolute value), | |, |δ|, |γ | ∈ [0.1, 0.15] as moderate elastic anisotropy and | |, |δ|, |γ | ∈ [0.15, +∞] as strong elastic anisotropy. The choice of these intervals is based on the divergence between the exact and approximate (by using Thomsen parameters) equations for the P-wave phase velocities in cracked media. Figure 8 shows the Thomsen-type anisotropy parameters in the Y -Z and X-Z planes (formulas are given in Appendix B). In the high-frequency limit, all anisotropy parameters are the same for both models with connected and disconnected cracks. Furthermore, for the model with disconnected cracks, all anisotropy parameters are frequency independent, because the stiffness tensor is frequency independent. For the model with connected cracks, several anisotropy parameters are frequency dependent due to squirt flow.
In the Y -Z plane, parameters Y Z and γ Y Z are zero for both models. The parameter δ Y Z is frequency dependent and controls the shape of the P-wave phase velocity between 0 and 90 • . In the high-frequency limit, δ Y Z exhibits the maximum negative value which corresponds to strong elastic anisotropy. As frequency decreases, δ Y Z also decreases reaching a zero value around 10 4 Hz, and then δ Y Z increases reaching its positive maximum at low frequencies, which corresponds to weak elastic anisotropy; the positive maximum is approximately one-third of the absolute value of its negative maximum. It is interesting that δ Y Z changes sign from negative to positive, which is indeed observed in the P-wave velocity behavior (Fig. 5a) as P-wave velocity changes polarity with frequency. This was also observed by Barbosa et al. (2017) in the framework of Biot's theory. This polarity change has a fully mechanical nature. In the high-frequency limit, cracks behave as hydraulically isolated and fluid highly stiffens the normal compliance of the cracks (not tangential). As frequency decreases, fluid started to flow from more compliant to stiffer cracks as a response to the applied displacement boundary condition. δ (yz) = 0 corresponds to zero anisotropy; the numerator of δ Y Z is [c 23 (ω) + c 44 (ω)] 2 − [c 33 (ω) − c 44 (ω)] 2 (see Ap-pendix B). Therefore, for zero anisotropy, c 23 (ω) + c 44 (ω) must be equal to c 33 (ω) − c 44 (ω). The function c 44 (ω) is constant across the whole frequency range, c 23 (ω) is strictly decreasing with frequency and c 33 (ω) is strictly increasing with frequency ( Fig. 4). At a certain frequency (here it is at ∼ 10 4 Hz), the c 33 and c 23 components are in such a combination that c 23 (10 4 ) + c 44 (10 4 ) ≈ c 33 (10 4 ) − c 44 (10 4 ), so δ Y Z = 0, and the P-wave velocity in the Y -Z plane behaves as in a fully isotropic media.
In the X-Z plane, XZ and δ XZ are frequency dependent in the model with connected cracks. XZ exhibits moderate elastic anisotropy at low frequencies, while δ XZ exhibits moderate elastic anisotropy at high frequencies. Other parameters are frequency independent and exhibit certain nonzero values from weak to moderate elastic anisotropy.

The universal elastic anisotropy index
The universal elastic anisotropy index A U (Ranganathan and Ostoja-Starzewski, 2008) is widely used to measure the anisotropy strength in crystallography, engineering and materials science. This parameter is designed to evaluate the anisotropy strength of crystals having any elastic symmetry class (Ranganathan and Ostoja-Starzewski, 2008). Since A U is a scalar, it gives a simple and fast identification of the overall anisotropy strength of a model.    shear and the bulk contributions simultaneously. In analogy to the universal elastic anisotropy index, two other parameters are adopted which define the anisotropy strength in bulk A bulk and in shear A shear . As far as we are concerned, these parameters have not been widely used in earth sciences: only a few studies were found. Almqvist and Mainprice (2017) applied the universal elastic anisotropy index and two similar parameters for bulk and shear to study seismic properties and anisotropy of the continental crust. Kube and De Jong (2016), Duffy (2018) and Vieira et al. (2019) applied A U to quantify the elastic anisotropy of polycrystals. A brief review of these anisotropic measures and all necessary equations for their calculation are provided in Appendix C. Figure 9 shows the universal elastic anisotropy index A U and the anisotropy measures in bulk A bulk (ω) and shear A shear . For the model with disconnected cracks, A U is constant and frequency independent (Fig. 9, black line). Because A U has a certain small value (about 0.058), the model with disconnected cracks exhibits a certain small anisotropy. For the model with connected cracks, A U in the high-frequency limit is almost the same as for the model with disconnected cracks (Fig. 9, red line). (The nature of the discrepancy is related to the region containing the crack intersection.) For the model with connected cracks, the overall anisotropy slightly decreases towards lower frequencies until 10 4.3 Hz, reaching its minimum of 0.048 (Fig. 9, red line). This local minimum indeed corresponds to the c 23 attenuation peak (Fig. 4c). Then, still towards lower frequencies, A U (ω) increases reaching its maximum of 0.083 at frequencies below 10 3 Hz (Fig. 9, red line). Thus, the overall anisotropy of the model mainly increases due to squirt flow between the cracks, so the crack connectivity increases the overall anisotropy of the model towards low frequencies.
The anisotropy measure in bulk A bulk is constant and frequency independent for the models with connected and disconnected cracks (Fig. 9b). It means that fluid flow does not affect bulk properties of the model or the anisotropy strength in bulk. On the other hand, the anisotropy measure in shear A shear (ω) basically reproduces the behavior of the universal elastic anisotropy index measure A U . Therefore, one can conclude that the fluid flow changes anisotropy in shear mode but not in bulk mode.

Elastic anisotropy
Thomsen-type anisotropic parameters provide a very detailed description of the velocity anisotropy in different planes. Most importantly, only a limited number of the stiffness tensor coefficients are needed to calculate , δ and γ in each plane. Thus, Thomsen parameters can be used to quantify the medium anisotropy using seismic data. On the other hand, when all components of the stiffness tensor are known and the model's symmetry is low, it is difficult to analyze the overall anisotropy due to a large number of Thomsen parameters. For example, if the model exhibits orthorhombic symmetry, one should analyze nine Thomsen anisotropic parameters (three in each plane). Due to a large number of Thomsen parameters in this study, it is difficult to evaluate whether the overall medium's anisotropy is increasing or decreasing with frequency and how far the current model is from the closest isotropy. Thus, in addition to (or instead of) the Thomsentype anisotropic parameters, the universal elastic anisotropy index can be used. The universal elastic anisotropy index A U and the related measures in bulk, A bulk , and shear, A shear , provide the overall description of the anisotropy strength regardless of the model's complexity. The calculation of these parameters is as simple as the calculation of the Thomsen parameters. An obvious disadvantage of the universal elastic anisotropy index (and related measures) is that it requires knowledge of the full stiffness tensor. Thus, this anisotropic measure can be useful to evaluate results of numerical simulations of laboratory experiments and for measuring the anisotropy of single crystals. The analysis of two sets of anisotropic measures shows that (i) the overall anisotropy of the model with connected cracks (Fig. 1) mainly increases due to squirt flow towards low frequencies with a slight local decrease at intermediate frequencies (Fig. 9a); (ii) in the Y -Z plane, the magnitude of the "delta" anisotropy parameter decreases, reaches zero, and then increases again (reaching approximately one-third of its high-frequency value) towards low frequencies (Fig. 8a, blue curve); and (iii) in the X-Z plane, the "delta" anisotropy parameter decreases towards low frequencies (Fig. 8b, blue curve), while the "epsilon" anisotropy parameter increases (Fig. 8b, green curve).

Comparison with previous works
In this study, we numerically solve a coupled fluid-solid deformation problem at the pore scale. If we consider the mesoscopic scale scenario and use Biot's (1941) equations, the fluid flow effects on the effective moduli are equivalent to that of the coupled elastic Stokes equations (as in the present study), as it was shown by Quintal et al. (2016). The frequency-dependent anisotropy due to fluid flow at the mesoscopic scale for orthogonal fracture sets with different degrees of connectivity was numerically studied by Rubino et al. (2017), but their study was limited to two dimensions. The main conclusion of Rubino et al. (2017) is that the anisotropy decreases with fracture connectivity in the seismic frequency band due to fluid flow between connected fractures. The Y -Z plane in the present 3D numerical model is reasonably equivalent to the 2D numerical model of Rubino et al. (2017) as well as the physical mechanism under consideration. The results of Rubino et al. (2017) are reflected in Fig. 5 of this study. However, a more in-depth analysis shows that the anisotropy in the Y -Z plane decreases, reaches zero and then increases again towards low frequencies due to squirt flow (Fig. 5a, green, yellow and blue curves, respectively). Moreover, our present study shows that the overall anisotropy of the model with cracks of finite length actually increases due to fluid flow between interconnected cracks (Figs. 8,9). This conclusion is not universal and is valid only for a specific set of model parameters. Barbosa et al. (2017) performed a more detailed study of seismic anisotropy for a similar fracture geometry in two dimensions, as in the study of Rubino et al. (2017), specifying that the decrease in anisotropy is described by the anisotropy parameter δ while is zero. Furthermore, they observed a polarity change of the P-wave phase velocity behavior with frequency. In the present study, the "delta" anisotropy parameter in the Y -Z plane is more pronounced in the low-frequency limit (Fig. 5a, blue curve) compared to the work of Barbosa et al. (2017) due to different material properties and the threedimensional nature of the present model configuration.
In summary, fluid flow effects on seismic anisotropy are nonlinear with a possible increase or decrease in the elastic anisotropy at different frequencies. These two extreme cases, the maximum negative and the maximum positive δ parameter (and hence P-wave velocity) in the Y -Z plane, correspond to the relaxed and unrelaxed states. In other words, seismic anisotropy may behave completely different in different scenarios; therefore, more studies should be performed, especially with the sensitivity analysis of model parameters.

A qualitative comparison against analytical models
Numerical simulations are useful but analytical models are especially attractive since they help us to better understand the physics and do not require sophisticated numerical simulations. The limitations of the analytical solutions are restricted to simple pore space geometry, and some assumptions related to physics are needed to derive the closed form analytical formulas. Such a comparison of the numerical results against an analytical solution has been performed by Alkhimenkov et al. (2020) for a different pore space geometry. Unfortunately, there is no analytical solution for the present study considering a periodic distribution of intersecting cracks in three dimensions. But the qualitative com-parison of the low-and high-frequency limits (which correspond to relaxed and unrelaxed states) is possible (Mavko and Jizba, 1991). Several analytical studies show that the anisotropy (described by Thompson's parameters) is, in general, more pronounced at high frequencies than at low frequencies Sarout, 2009, 2011). In the relaxed state, one can calculate the effective dry elastic moduli and use Gassmann's equations to obtain the effective moduli of the saturated medium. In the unrelaxed state, one can calculate the effective elastic moduli by restricting fluid flow (by using zero displacement boundary conditions in the crack intersections). The low-and high-frequency limits for elastic moduli have been calculated using these semi-analytical approaches, and numerical results have been reproduced.

Conclusions
We have numerically calculated the frequency-dependent elastic moduli of a fluid-saturated porous medium caused by squirt flow. The considered 3D numerical models consist of two perpendicular connected or disconnected cracks embedded in a solid grain material. Cracks are represented by very flat cylinders filled with a fluid. Grains are described as a linear isotropic elastic material, while the fluid phase is described by the quasistatic linearized compressible Navier-Stokes momentum equation.
We showed that seismic velocities are azimuth, angle and frequency dependent due to squirt flow between connected cracks. The resulting elastic frequency-dependent anisotropy was analyzed by using the Thomsen-type anisotropic parameters and the universal elastic anisotropy index. The latter is a scalar parameter which can be used to analyze the overall anisotropy of the model and its divergence from the closest isotropy. We showed that the seismic anisotropy may locally decrease as well as increase due to squirt flow in one specific plane. However, the overall anisotropy of the model mainly increases due to squirt flow between the cracks towards low frequencies. In the model with disconnected cracks, no fluid flow occurs, and thus the effective properties of the model correspond to the elastic limit. The elastic limit is equivalent to the high-frequency limit for the model with connected cracks. Seismic velocities are only azimuth and angle dependent as for a fully elastic material, and they are independent of frequency.
In summary, squirt flow does affect effective mechanical properties of cracked rocks and thus seismic velocity anisotropy. Given that seismic anisotropy variations with frequency are very sensitive to the pore space geometry and material properties, it is difficult to make a general prediction. According to our study, the effective frequency-dependent response of a cracked medium is different in different planes. The local response (in a certain plane) is controlled by crack orientation, which is the key parameter. The magnitude of the frequency-dependent response (i.e., the dispersion and attenuation) is controlled by crack compliances, crack porosity and their fluid content. (Dry or liquid-saturation conditions will cause completely different behavior.) Most importantly, crack porosity is a very important parameter in fluidsaturated rocks (contrary to dry rocks) since it defines the volume of fluid which may flow due to wave propagation, causing wave attenuation and dispersion.
The mixed test for the c 13 component can be derived from the anisotropic stress-strain relation (Hooke's law) (similarly to the c 13 component in a VTI medium; Alkhimenkov et al., 2020); xyL is set to u zz = u; u xx and u yy are free; yzL is set to u xx = u; and u zz and u yy are free, where u = 10 −6 . In the other four planes, the normal component of the displacement is set to zero and other components are free. The the stress-strain relation for the σ 1 stress component is Using the abovementioned boundary conditions and setting 2 = 0, Eq. (A1) becomes Equation (A2) can be solved for the c 13 component; the solution is Equation (A3) is used to calculate the c 13 component (c 11 is taken from the direct tests). The mixed test for the c 23 (in this numerical model, the c 23 is dispersive) component again can be derived from the anisotropic stress-strain relation (Hooke's law) (similarly to the previous test). xyL is set to u zz = u; u xx , u yy are free; xzL is set to u yy = u; and u zz and u xx are free.
In the other four planes, the normal component of the displacement is set to zero, other components are free. Then, using the following equation the c 23 component is calculated (c 22 is taken from the direct test). Equations (A3)-(A4) are found from Hooke's law considering nonzero strains in the x and z (in y and z) directions and then solving a system of two equations analytically.

Appendix B: Thomsen-type anisotropic parameters
Thomsen-type anisotropic parameters (Thomsen, 1986) are widely used in the applied geophysics community. Thomsen weak anisotropy parameters were originally developed for vertical transverse isotropic media (Thomsen, 1986). A natural extension of these parameters to orthorhombic media was suggested by Tsvankin (1997) and Bakulin et al. (2000a). These parameters correspond to the anisotropy of the compression and shear waves in orthorhombic media in different Cartesian propagation planes. In the Y -Z plane, Thomsentype anisotropic parameters are Assuming that one deals with an anisotropic frequencydependent effective fourth-rank stiffness tensor C (might be frequency dependent C ⇒ C(ω)), a compliance tensor is defined as S(ω) = C(ω) −1 . Then, for each frequency, the effective single orientation fourth-rank stiffness and compliance tensors are uniformly distributed and the isotropic stiffness and compliance tensors are calculated. Averaging the single orientation stiffness tensor belongs to the Voigt assumption, which is the theoretical maximum value of the isotropic elastic moduli. On the other hand, averaging the single orientation compliance tensor belongs to the Reuss assumption which is the theoretical minimum value of the isotropic elastic moduli. The resulting isotropic tensors can be expressed in terms of the spherical and deviatoric parts corresponding to bulk K and shear moduli µ: and where the superscripts "V" and "R" correspond to Voigt and Reuss estimates, respectively. J and D are the spherical (volumetric) and deviatoric parts of the symmetric unit fourthorder tensor. The double contraction of the scalar product (quadruple contraction) of Eqs. (C1) and (C2) gives If the effective stiffness tensor is isotropic, then C V (ω) = S R (ω) −1 and K V /K R = µ V /µ R = 1. Therefore, when the effective stiffness tensor is isotropic, the value of Eq. C3 equals to 6 and this value increases when the effective stiffness tensor becomes anisotropic. Thus, the universal elastic anisotropy index measure A U is defined as (Ranganathan and Ostoja-Starzewski, 2008): In geophysics, the separation of the elastic anisotropy measure in bulk and shear modes is necessary because rocks might exhibit different frequency dependencies due to bulk and shear deformations. Therefore, in analogy to the universal elastic anisotropy index measure A U , the anisotropy measures in bulk A bulk (ω) and shear A shear (ω) can be defined as and These two parameters, A bulk (ω) and A shear (ω), obey the same interpretation as the universal elastic anisotropy index measure: A bulk = 0 (A shear = 0) corresponds to zero bulk (shear) anisotropy of the model, while the discrepancy of A bulk (A shear ) from zero anisotropy defines the anisotropy strength in bulk (shear).
The Voigt and Reuss estimates (K V , K R , µ V and µ R ) can be calculated via simple algebraic formulas for different symmetry classes which can be found elsewhere, e.g., in Ravindran et al. (1998) for orthorhombic symmetry, in Feng et al. (2012) for tetragonal symmetry and in Duffy (2018) for cubic symmetry. Thus, for orthorhombic symmetry (lowest possible symmetry created by two perpendicular sets of cracks embedded into an isotropic material), the Voigt and Reuss bulk moduli are (written for the components of stiffness (compliance) matrices c ij (s ij ), in Voigt notation) Equations (C7)-(C10) are valid for orthorhombic symmetry and for higher symmetries: tetragonal, transverse isotropy and cubic. Thus, for evaluating the universal elastic anisotropy index A U and the anisotropy measures in bulk A bulk (ω) and shear A shear (ω), one can use Eqs. (C7)-(C10) to calculate the Voigt and Reuss estimates (K V , K R , µ V and µ R ) and then calculate A U using Eq. (C4) and A bulk (ω) and A shear (ω) using Eqs. (C5) and (C6), respectively.