Impact of permeability anisotropy misalignment on flow rates predicted by hydrogeological models

The notion of permeability is critical to compute underground fluid flow. In most cases rock permeability is anisotropic, due to physical processes including gravitational compaction, which often results in the principal permeability directions being approximately horizontal and vertical in undeformed rocks. However, rocks often are tilted and/or deformed over time, therefore permeability orientation varies. Anisotropic permeability with varying orientation is hard to quantify in three-dimensional (3D) models and is therefore sometimes approximated, for convenience, by setting the principal permeability directions to horizontal and vertical, and assuming that corresponding errors in fluid flow might be negligible when the change in orientation is minimal. This study shows how minor misalignment of the permeability tensor can lead to large errors in fluid flow magnitude and corresponding transport times for strongly anisotropic rocks. It also provides a method to set anisotropic permeability orientation appropriately in geometrically complex 3D models using implicit 3D geological modelling. The misalignment is particularly costly when fluid flow is localised in thin channels, where a misalignment of just 5° leads to errors of two orders of magnitude for anisotropy ratios (between the largest and smallest principal values of the permeability tensor) of 104. It is therefore recommended to set anisotropic permeability accurately, using longitudinal and transverse components along with their respective orientations, rather than horizontal and vertical components. This approach will become increasingly important as 3D models gain realism in their representation of complex geometries.


Introduction
Permeability is a key controlling factor for fluid flow in the subsurface (Bjørlykke 1993).Permeability anisotropy exists, especially in sedimentary rocks, due to various processes including gravitational compaction and grain alignment.These processes result in permeability typically being larger parallel to bedding than perpendicular to it (Bear 1972).
Permeability anisotropy in sedimentary rocks is characterised by the ratio R = k ∥ /k ⊥ , where k ∥ and k ⊥ denote the permeability components parallel (longitudinal) and perpendicular (transverse) to bedding.Gravitational compaction of an isotropic pore network can result in R values around one order of magnitude (Scholes et al. 2006); however, real values vary by more than four orders of magnitude for rocks with the same porosity, mean pore throat size and clay mineral content (Armitage et al. 2011).Tilting and folding of sedimentary rocks lead to spatial variations in the orientation of k ∥ and k ⊥ that can have a significant impact on fluid flow.
Capturing the direction of permeability anisotropy in nonplanar sedimentary rocks is critical for accurate prediction of fluid flow.This requires a fluid flow simulator that includes full tensor representation of permeability.While this functionality exists in several fluid flow simulators-e.g.MOOSE (Wilkins et al. 2020); SUTRA (Provost and Voss 2019); MODFLOW 6 (Langevin et al. 2017); FEMWATER (Lin et al. 1997)-it is difficult to use without an automated method for populating the permeability tensor to reflect the varying orientation of strata throughout the modelled domain.Furthermore, simulators that do not provide a full tensor representation of permeability such as the earlier versions of MODFLOW, are still widely used (Bardot et al. 2022).Consequently, the principal permeability directions are commonly approximated as being horizontal and vertical in simulations of subsurface fluid flow, or the coordinate axes are aligned with the dominant principal permeability directions, which still does not capture local variations in permeability orientation.To understand the errors that may be introduced by this approximation, consider the case of fluid flowing along a slanted rock layer dipping at angle α with anisotropy ratio R, and with k ∥ and k ⊥ parallel and perpendicular to the layer (Fig. 1a).A fluid pressure gradient between the ends of the layer drives fluid flow parallel to the layer; however, if the principal permeability directions are assumed to be horizontal and vertical the fluid flow direction will be at an angle to the layer (Fig. 1b), and the resulting error in the fluid transport time can be nonnegligible depending on R and the geometry of the layer.To understand this, consider the case where R → ∞.In such a case, a slight misalignment between the true and modelled permeability directions may be enough to block completely the flow, assuming the layers above and below are impermeable and the layer is long enough that there is no horizontal path from one end to the other.
This simple example demonstrates the importance of considering a proper orientation for the permeability tensor when the bedding is not perfectly horizontal, including where the bedding orientation varies spatially (as in folded sedimentary rocks) and/or between different geological units.This is particularly important when transport times and flow rates are at the core of the problem, as in the fields of contaminant transport, waste storage, hydrocarbon/ geothermal reservoir engineering, hydrogeology, and the formation of mineral deposits.For example, Yager et al. (2009) showed that different representations of permeability anisotropy resulted in significantly different groundwater ages and fluid flow pathways in models of a folded aquifer system in the Shenandoah Valley, USA.Similarly, Bardot et al. (2022) demonstrated significant differences in fluid flow directions predicted by a hydrogeological model of the Perth Basin (Western Australia) when assuming horizontal/vertical permeability directions rather than longitudinal/ transverse directions aligned with the dipping aquifer units.In a different application, Poulet et al. (2022) demonstrated the importance of accurate representation of permeability anisotropy for understanding the distribution of supergene iron ore deposits in folded sedimentary rocks.
The purpose of this paper is two-fold: (1) to quantify the errors in fluid flow rate that arise due to incorrect representation of the permeability tensor, as a function of R and the geometry of a permeable layer; and (2) to present a method for populating the permeability tensor appropriately in 3D models using implicit 3D geological modelling, such that the errors previously described can be avoided.

Methods
This section presents a methodology to set the permeability tensor appropriately in a geological model, i.e. following the geological structures affecting fluid flow, such as sedimentary bedding or faults.The workflow (Fig. 2)

a b
Fig. 1 Conceptual comparison of fluid flow in slanted permeable layers, at an angle α from the horizontal, when the permeability tensor is a aligned with the layers or b not aligned with the layers.All other parameters are identical, including the driving force imposing a gradient of pressure from left to right.The Darcy flow q is aligned with the layer and will therefore remain constant for a given k ∥ regardless of the value of R (a).The misalignment between the layer and the permeability tensor forces the Darcy flow q at a shallower angle (b).For a thick (bottom) layer, there might still be a channel available at that angle between the two sides of the model, yet reduced in size (indicated by pink shading).For a thin (top) layer, only the sides allow the flow to follow this preferred direction over a short distance.For the main part of the thin layer, the flow is forced in the direction of the layer and will therefore have a reduced magnitude as R increases extends to three dimensions of the process described in (Poulet et al. 2022).The description is focused here on following the strata in a sedimentary basin, although the method is applicable to other volumetric geological features such as fault zones or igneous bodies.
The process is based on digitised surfaces selected by the modeller (Fig. 2a).They should include the boundaries between the geological layers of interest, but can also include extra surfaces within the layers to reflect, for instance, the orientation of bedding within major units.Such surfaces can be exported from traditional geological modelling packages or generated parametrically for synthetic models.Those surfaces are used to generate both a potential field to align the permeability appropriately, and a mesh that will be used to run numerical fluid flow simulations.Among the various approaches available for the potential field, implicit modelling is particularly well adapted (Wellmann and Caumon 2018).This study uses the open-source implicit modelling package LoopStructural (Grose et al. 2021) to generate that potential field (Fig. 2b), in which all points on each of the provided surfaces have a common value, and values in between the surfaces are obtained by interpolation, noting that the code handles discontinuities (e.g.faults, unconformities).The gradient of the potential field defines the direction of transverse permeability (k ⊥ ), with longitudinal permeability (k ∥ ) assumed to be equal in the orthogonal directions.Note that absolute values of the potential field are not relevant since it is only the orientation of the isosurfaces of this field that are required to align the permeability.
While implicit modelling packages could theoretically generate the corresponding mesh for fluid flow simulations, their mesh generation functionalities are usually not as welldeveloped as dedicated meshing packages; therefore, a separate meshing package, gmsh (Geuzaine and Remacle 2009), is used to create a mesh matching the same provided surfaces (Fig. 2c).Any discrepancy stemming from the use of two separate methods for the mesh and potential field can be controlled through the resolution of the digitised input points.
The computation of the permeability tensor Κ at any point within the model is computed from the user-defined values of longitudinal (k ∥ ) and transverse (k ⊥ ) permeability, following the relationship where t is the unit normal vector of the potential field and I the identity tensor.In the selected example, where the anisotropy follows the stratigraphy, values of k ∥ and k ⊥ are set separately for each geological unit.This permeability computation (Fig. 2d) is implemented in the Redback opensource simulator (Poulet and Veveakis 2016), which is based on the finite element MOOSE framework (Permann et al. 2020).The resulting permeability field and corresponding mesh can then be used to simulate geological fluid flow with inert tracer transport (Fig. 2e) for instance, here using the open-source PorousFlow module (Wilkins et al. 2020) of MOOSE.
This workflow to populate the anisotropic permeability field using implicit geological modelling is not required for all fluid flow simulators.For example, in some simulators (e.g.ECLIPSE TM ; Schlumberger, 2014) the principal permeability directions can be parallel to the local coordinate axes of the mesh, which may follow geological layering.
(1) The methodology described in this article is useful for fluid flow simulators that do not provide this functionality, providing greater flexibility to simulate fluid flow on meshes that do not necessarily follow geological layering.
The impact of varying the direction of permeability anisotropy is illustrated by simulating single-phase fluid flow in fully saturated porous rock by solving the continuity equation where ϕ denotes the porosity, ρ f the fluid density (kg m −3 ), and v f the Darcy fluid flow velocity (m s −1 ), expressed as with k the permeability tensor (m 2 ), μ the fluid viscosity (Pa s), P the pore fluid pressure (Pa), and g (m s −2 ) the gravitational body force.Mass (tracer) transport is not simulated in the examples presented in this study.However, the implications of the results for tracer transport can be readily inferred-for example, a reduction in fluid velocity due to misalignment of the permeability would result in a tracer travelling a shorter distance in a given time, or equivalently, the transport time over a given distance would be longer.

Application to folded strata
The approach is illustrated in a realistic geometry including faults and folded sedimentary rocks.The example is based on the Sheep Mountain Anticline (Fig. 3a), where fluid flow is simulated from the base of the model through a network of (2) faults more permeable than the host rock, to reach a permeable limestone layer in between two shale layers with lower permeability.Three permeability scenarios were tested: (1) isotropic permeability; (2) anisotropic permeability with principal directions being horizontal and vertical; and (3) anisotropic permeability with principal directions following the orientation of the geological units, using the implicit geological modelling method described in the preceding to set the appropriate permeability orientations.The flow pattern within the faults is not of interest in this case and the fault permeability is therefore considered to be isotropic in all scenarios; anisotropy in faults is considered in the discussion section below.it reaches the limestone layer, where the fate of the fluid depends on the permeability scenario.Isotropic permeability (Fig. 3b) results in pervasive flow through the limestone and overlying shale layer, allowing fluid to escape through the top boundary.Anisotropic permeability with horizontal/ vertical orientation (Fig. 3c) results in a dominantly horizontal flow direction in the limestone and overlying shale, again allowing fluid to escape through the shale to reach the top boundary.In the most realistic case, where the principal permeability directions follow the geological layers (Fig. 3d), the flow remains largely confined within the limestone unit.The patterns of Darcy velocity magnitude in the limestone highlight the ability to capture patterns of fluid flow around fault splays within a geological unit through the implicit modelling of the potential function matching the complex geometry of the folded sedimentary rocks.In all cases, there is a thin region of high Darcy velocity at the base of the model-this is an artefact of the elevated fluid pressure boundary condition there, and is not representative of the true fluid flow regime in the low permeability basement.Note that the values of longitudinal and transverse permeabilities used in this example are generic values for the relevant rock types, not values specific to the modelled location.

Effect of permeability misalignment in a dipping layer
To verify the underpinning argument for this study, the effect of misalignment of the permeability is quantified by simulating some scenarios based on Fig. 1, for the case where the permeability tensor is not correctly aligned with the strata.A generic two-dimensional (2D) dimensionless model is used to represent a single permeable layer (with anisotropic permeability, k h /k v = R), oriented at an angle α from the horizontal in a less permeable host rock with isotropic permeability value k h /10 7 (Fig. 4).The permeability tensor is kept aligned with the horizontal and vertical axes of the model.For comparison purposes, the dimensionless thickness h of the central layer is fixed arbitrarily to 0.01, the dipping angle α to 5 °, and various values of R and the dimensionless length L are considered.For a short length (L = 0.05), the geometry is such that the inlet and outlet are connected directly by a fluid pathway in the preferred flow direction, like in the case of the thick layer of Fig. 1b.Longer model scenarios (L = 0.2, L = 1) are also used, for which the inlet and outlet are no longer directly connected for some values of anisotropy ratio R, like in the case of the thin layer of Fig. 1b.Structured meshes comprising quadrilateral cells are used in all three scenarios, with 40 cells vertically in the central layer for L = 1 and 100 cells for the other two cases (L = 0.2 and L = 0.05), as longer meshes require more elements to reach the same numerical accuracy (mesh dependency analyses not shown).The corresponding meshes contain 60,000 cells (for L = 0.05) or 240,000 cells (for L = 0.2 and L = 1), with the horizontal resolution increasing towards the left and right boundaries (Shishkin mesh, Kopteva and O'Riordan 2010) to capture sharp changes near those boundaries (boundary layers).The meshes are created such that the edges of the central layer coincide with cell edges, to avoid stair-step effects associated with orthogonal meshes.A gradient of pore pressure is imposed between the left and right vertical boundaries of the central layer with Dirichlet (i.e.fixed pore pressure) boundary conditions, with no flux on the top and bottom boundaries.
An isothermal pore pressure diffusion problem is then solved using the open source PorousFlow module (Wilkins et al. 2020) of the finite element MOOSE simulator (Permann et al. 2020).For each value of L, the magnitude of the fluid flow within the central layer is investigated for various values of anisotropy ratio (R ∈ {1, 10, 10 2 , 10 3 , 10 4 }), by tracking the fluid flow magnitude at the centre of the model (point C in Fig. 4).The resulting distribution of fluid flux magnitude is illustrated in Figs. 5, 6 and 7.The fluid flow direction and magnitude predicted for the isotropic case (R = 1) is the same as it would be with R > 1 if the permeability tensor were oriented "correctly" (i.e. with principal permeability directions being parallel and perpendicular to the layer, rather than horizontal and vertical); hence, the isotropic case is used to quantify the error in fluid flow magnitude that arises from misalignment of the permeability tensor when R > 1. Note, however, that the equivalence of the isotropic case to the anisotropic case with correctly aligned permeability is a special feature of this simple scenario; in more complex 3D geometries where the orientation of layers varies spatially (e.g. in folded layers), and/or where there is movement of fluid between layers, the isotropic and correctly oriented anisotropic cases would produce different results.
The results, reported in Table 1, are normalised by the fluid flow magnitude in the isotropic case (R = 1) to assist  7) rather than following the dip of the layer.When the inlet and outlet are not connected along the preferential flow direction-L = 0.2 (Fig. 5), L = 1 (Fig. 6)-the error in flow magnitude in the central layer is already ∼6% for a moderate anisotropy ratio R = 10.This error becomes overwhelming in cases of large anisotropy, with fluid flow reduced by two orders of magnitude when R = 10 4 .Even in the case L = 0.05, where a conduit connects the inlet and outlet following the preferred flow direction, while the amplitude of the flow in that conduit is not affected, the thickness of that conduit is less than that of the layer (Figs. 5 and 7), thus reducing the overall amount of whatever quantity that would be advected by fluid flow.

Discussion
Neglecting anisotropy orientation when modelling fluid flow in sedimentary rocks can lead to considerable fluid flow magnitude discrepancies (Table 1) even at low strata angles, with the advection of any other property being automatically affected by that flow magnitude.This can lead to critical transport timing errors for geological applications like contaminant transport or nuclear waste disposal for instance, as highlighted in Fig. 6 by the two orders of fluid flow magnitude error for high anisotropy ratio, translating directly to the same error in advection velocity.In particular, for nuclear waste disposal, discrepancies in the computed flow fields will result in significant errors in estimated radionuclide transport times/ distances, as the safety assessment of radioactive waste repositories is typically conducted for a time span of 1 million years (Metcalfe et al. 2008).Slight undulations or imperfect parallelism of layers could result in significant errors even where the layers appear to be largely planar and parallel on a large scale, if the anisotropy direction is not accurately represented in the model.It is therefore preferable to avoid the use of horizontal and vertical permeabilities (with ratio k h /k v ) and use instead the notions of longitudinal and transverse permeabilities (with ratio k ∥ /k ⊥ ) with their corresponding orientations, even when the two sets might seem almost identical.This is particularly important when the anisotropy ratio is large.There is a computational cost associated with the approach proposed, mostly in terms of memory required to store the permeability tensor field across the model, but this cost should be assessed in light of the improved accuracy of the corresponding results.This applies both for simulators requiring a preprocessing step to populate the anisotropic permeability field (e.g. the implicit geological modelling   4 with α = 5 °, L = 1 and h = 0.01, using the same colour scheme as in Fig. 4. The results match those obtained for L = 0.2 (Fig. 5), only with the geometrical effects at the side boundaries affecting less the overall behaviour (away from those boundaries) methodology presented in the section "Methods"), as well as simulators that provide the functionality to align permeability with the geological layers, e.g.ECLIPSE TM (Schlumberger, 2014).
An alternative approach that avoids the requirement to specify anisotropy is to divide each layer into many small sublayers, each sublayer having isotropic permeability, but with contrasting values of permeability between adjacent sublayers.This could represent alternating layers of varying grain size within a sedimentary unit, for example.This explicit representation of sublayers will result in the overall effective permeability of the layer being anisotropic; however, it comes at significant computational cost due to the large number of mesh elements required to represent the thin sublayers, and fails to capture anisotropy within the sublayers, which would be expected in clay-rich sublayers, for example, due to grain alignment.Thus, the proposed approach is a useful method for homogenising and upscaling the permeability of a rock unit in which anisotropy exists on various scales.The upscaled anisotropic permeability can account for both the layering of different sediment types with contrasting permeability (which may in itself be anisotropic) and the connectivity across those sublayers due to gaps, which leads to smaller anisotropy ratios at larger scales.The use of potential fields to capture the geometry of complex geological structures can also account for the superposition of fracture-related permeability that occur in folded strata, as it was shown that the fold geometry strongly dictates the arising fracture pattern (Ramsay and Huber 1987).
Note that for cases where thin permeable layers surrounded by lower permeability rocks are explicitly represented in the model, it may be acceptable to use isotropic permeability equal to the longitudinal permeability for those thin layers if the flow direction is expected to be predominantly along the layers; indeed, this would be preferable to assigning anisotropic permeability with incorrect orientation.Similarly, if flow is expected to be predominantly across the layer, it may be appropriate to assign isotropic permeability equal to the transverse permeability.However, anisotropic permeability would be required if it were important to resolve the exact fluid flow pathway through the layer, e.g. if a tracer propagation front in that layer needs to be resolved accurately, especially if there are significant variations in orientation along the layer.It is only in very simple scenarios, such as that represented in Fig. 3, that the isotropic case produces identical fluid pathways to the anisotropic case with correctly aligned anisotropy.For lower values of anisotropy, here shown with R = 10, the geometrical effects around the side boundaries visibly affect the flow and distort the streamlines.To highlight the variations responsible for that effect, the colorbar range for the case R = 10 is reduced to two orders of magnitude Table 1 Normalised fluid flow magnitude, as a percentage of the isotropic case, at the centre of the model (point C on Fig. 4) for various values of model length L and anisotropy ratio R, at fixed angle α = 5° and layer thickness h = 0.01.For L = 0.05, the geometrical effects near the side boundaries are quite pronounced for lower values of R and streamlines are not parallel to the central layer (Fig. 7).This effect explains why the normalised flow magnitudes are not exactly 100% for R = 10 and R = 10 2 for the choice of Dirichlet boundary conditions, despite the existence of a pathway between the inlet and outlet along the preferred flow direction in these cases

Anisotropy ratio
Model length Another alternative approach to generating a 3D model with correctly aligned anisotropic permeability uses stratigraphic forward modelling (SFM) to obtain a 3D distribution of porosity and sediment types by simulating deposition of sediments in a sedimentary basin.An upscaled anisotropic permeability tensor can then be obtained at each cell of a superposed mesh for fluid flow simulation, derived from the orientation, porosity and sediment composition of cells in the stratigraphic forward model, given some appropriate porositypermeability relationships (Sheldon et al. 2023).On the one hand, such an approach can provide a realistic anisotropic permeability distribution, depending on adequate calibration of the model, while on the other, SFM requires considerable skill, time, and data for calibration.It might also require an additional step to transform the generated mesh and permeability field if fluid flow needs to be simulated at a later time when the basin has been mechanically deformed.
Errors induced by the misalignment of anisotropic permeability are not limited to strata orientation, with other geological features playing equally important roles in controlling fluid flow.This is the case for instance with faults, which can be strongly anisotropic in terms of their permeability (e.g.Bense et al. 2013), with the direction of anisotropy following the nonplanar geometry of the fault zone.Capturing the appropriate local permeability tensor orientation in fault zones represented as volumetric features (i.e.features of finite width) in a 3D model can be addressed with the same methodology as described previously for sedimentary layers.Faults typically represent subvertical fluid pathways and/or barriers to cross-fault flow, and modelling their hydraulic behaviour accurately is important to complement the strata's typically subhorizontal fluid channelling ability in sedimentary basins.This can be particularly useful for understanding/predicting the formation of sediment-hosted mineral deposits for example, where the mineralising fluid moves between faults and sedimentary rocks, or for understanding groundwater flow paths where an aquifer is cut by faults.For these applications, the choice of approach for representing fault permeability depends on the scale and permeability characteristics of the faults and their orientation with respect to the dominant fluid flow direction-for instance, anisotropy may be associated with a preferred alignment of fractures in the fault zone, or by a low permeability fault gouge surrounded by a more permeable damage zone.In such cases, one could represent the fault zone as a volumetric feature of finite width, with appropriate anisotropic permeability to resolve fluid flow patterns within the fault zone if required; alternatively, if the thickness of the fault zone is very small relative to the overall size of the model, it may be represented as a lower-dimensional feature using the techniques described by Poulet et al. (2021) to capture the anisotropy imposed by each component of the fault zone (i.e.damage zone(s) and core).In other cases, anisotropy within a fault zone may not have a significant impact on the overall fluid flow regime, in which case the fault can be represented with isotropic permeability; for example, this may be the case when the dominant flow pathway involves fluid flowing up steep faults and then deviating out into subhorizontal sedimentary rocks (e.g.Sheldon et al. 2023).In such cases the transverse permeability of the fault is largely irrelevant; however, where the dominant fluid flow direction is at a high angle to faults (e.g.subhorizontal flow interacting with subvertical faults), it is important to capture the anisotropic permeability of the faults and their host rocks as this will have a strong influence on fluid pathways through the system (e.g.Bense and Person 2006).

Conclusion
In conclusion, this study demonstrates the importance of considering longitudinal and transverse permeability values and their corresponding orientations rather than assuming the principal permeability directions to be horizontal and vertical, which can lead to large errors in simulated fluid flow directions and magnitudes, and corresponding errors in transport times, as the anisotropy ratio grows.This is particularly important for applications like contaminant transport or nuclear waste disposal, for which quantifying travel times accurately is essential, but also plays a nonnegligible role for any fluid flow modelling application in the subsurface, including in the fields of mineral exploration, hydrocarbon reservoir modelling, geothermal energy, carbon sequestration, or hydrogeology.
The methodology described in this study addresses the difficulty of populating the anisotropic permeability field in 3D models of deformed geological units (e.g.faults and folds), by using implicit geological modelling to provide the orientation of the permeability tensor throughout the model.
The increased accuracy of simulated fluid flow fields that is achieved by accurate handling of anisotropic permeability justifies the added numerical cost in many cases, and can only become more important in future modelling studies with increasingly complex 3D geometries.

Fig. 2
Fig. 2 Workflow used to generate anisotropic permeability for fluid flow simulation.a A set of surfaces is used to generate both b a potential function over the whole volume and c a corresponding mesh.The gradient of the potential function defines the direction of transverse permeability, to produce a field of permeability ten- Pore pressure Dirichlet boundary conditions are imposed on all sides of the model, with a constant value on the top surface, equilibrated hydrostatic values on the sides after an initialisation stage, and excess pore pressure at the bottom to force the fluid upwards.These fluid flow boundary conditions are selected to drive upward fluid flow and evaluate the relative impact of each permeability scenario; they do not represent reality.A geothermal gradient of 25 °C is applied with a top temperature of 10 °C and the fluid properties vary with temperature and pore pressure following the 1997 implementation of the International Association for the Properties of Water and Steam (IAPWS 2014).Porosity and permeability values representative of the geological units are used throughout the model.The porosity is set to 1% (resp.10%) in the basement and shale (resp.limestone and fault) units.The permeability is set to 10 −17 m 2 in the basement, 10 −16 m 2 in the shale layers, 10 −14 m 2 in the limestone layer and 10 −13 m 2 within the faults.For the anisotropic scenarios, those values in the limestone and shale units become the horizontal or longitudinal components, while the vertical or transverse components are divided by a factor 10. Plots of the resulting fluid flow magnitude (Fig. 3b-d) indicate that flow is confined mainly to the faults until

Fig. 3
Fig. 3 Example inspired by the Sheep Mountain Anticline scenarios described in (Beaudoin et al. 2023).a Conceptual geological set up, where the faults and limestones have greater permeability values than other units (vertical exaggeration 1.3×).Darcy velocity magnitudes from indicative numerical simulations of fluid flow from the base of the model are shown for three scenarios, using b isotropic permeability, c anisotropic permeability with principal directions being horizontal and

Fig. 4
Fig. 4 Setup for 2D flow simulations on a single layer

Fig. 5
Fig. 5 Fluid flow magnitude for different anisotropy ratios R in two cases of the model from Fig. 4.Both have the same layer thickness h = 0.01, same angle α = 5 °, but different values of model length, L = 0.05 or L = 0.2.All Darcy velocity values are normalised by the maximum value per simulation, and the colorbar is in logarithmic scale between 10 −4 and 1.Those results reproduce the cases of the thick (L = 0.05) and thin layers (L = 0.2) in Fig. 1

Fig. 6
Fig.6Normalised fluid flow magnitude for different anisotropy ratios R in a case of the model from Fig.4with α = 5 °, L = 1 and h = 0.01, using the same colour scheme as in Fig.4.The results match those obtained for L = 0.2 (Fig.5), only with the geometrical effects at the side boundaries affecting less the overall behaviour (away from those boundaries)

Fig. 7
Fig. 7 Normalised fluid flow magnitude and streamlines in the central layer for R = 10 and R = 10 4 , in the case L = 0.05.For lower values of anisotropy, here shown with R = 10, the geometrical effects around the side boundaries visibly affect the flow and distort the streamlines.To highlight the variations responsible for that effect, the colorbar range for the case R = 10 is reduced to two orders of magnitude