Grid refinement for modeling multiphase flow in discretely fractured porous media

A study of the effects of grid discretization on the migration of DNAPL within a discrete-fracture network embedded in a porous rock matrix is presented. It is shown that an insufficiently fine discretization of the fracture elements can lead to an overprediction of the volume of DNAPL that continues to migrate vertically at the intersection of a vertical and horizontal fracture. Uniform discretization of elements at the scale of one centimetre (or less) accurately resolved the density and capillary pressure components of the head gradient in the DNAPL. An alternative, non-uniform method of discretization of elements within the discrete-fracture network is presented whereby only fracture elements immediately adjacent to fracture intersections are refined. To further limit the number of elements employed, the porous matrix elements adjacent to the fracture elements are not similarly refined. Results show this alternative method of discretization reduces the numerical error to an acceptable level, while allowing the simulation of field-scale DNAPL contamination problems. The results from two field-scale simulations of a DNAPL-contaminated carbonate bedrock site in Ontario, Canada are presented. These simulations compare different methods of grid discretization, and highlight the importance of grid refinement when simulating DNAPL migration problems in fractured porous media.     

Integrity of Clay Till Aquitards to DNAPL Migration: Assessment Using Current and Emerging Characterization Tools

Field investigations were carried out to determine the occurrence of tetrachloroethene (PCE) dense nonaqueous phase liquid (DNAPL), the source zone architecture and the aquitard integrity at a 30‐ to 50‐year old DNAPL release site. The DNAPL source zone is located in the clay till unit overlying a limestone aquifer. The DNAPL source zone architecture was investigated through a multiple‐lines‐of‐evidence approach using various characterization tools; the most favorable combination of tools for the DNAPL characterization was geophysical investigations, membrane interface probe, core subsampling with quantification of chlorinated solvents, hydrophobic dye test with Sudan IV, and Flexible Liner Underground Technologies (FLUTe) NAPL liners with activated carbon felt (FACT). While the occurrence of DNAPL was best determined by quantification of chlorinated solvents in soil samples supported by the hydrophobic dye tests (Sudan IV and NAPL FLUTe), the conceptual understanding of source zone architecture was greatly assisted by the indirect continuous characterization tools. Although mobile or high residual DNAPL (S _t > 1%) only occurred in 11% of the source zone samples (intact cores), they comprised 86% of the total PCE mass. The dataset, and associated data analysis, supported vertical migration of DNAPL through fractures in the upper part of the clay till, horizontal migration along high permeability features around the redox boundary in the clay till, and to some extent vertical migration through the fractures in the reduced part of the clay till aquitard to the underlying limestone aquifer. The aquitard integrity to DNAPL migration was found to be compromised at a thickness of reduced clay till of less than 2 m.     

The Effect of Modeled Recharge Distribution on Simulated Groundwater Availability and Capture

Simulating groundwater flow in basin‐fill aquifers of the semiarid southwestern United States commonly requires decisions about how to distribute aquifer recharge. Precipitation can recharge basin‐fill aquifers by direct infiltration and transport through faults and fractures in the high‐elevation areas, by flowing overland through high‐elevation areas to infiltrate at basin‐fill margins along mountain fronts, by flowing overland to infiltrate along ephemeral channels that often traverse basins in the area, or by some combination of these processes. The importance of accurately simulating recharge distributions is a current topic of discussion among hydrologists and water managers in the region, but no comparative study has been performed to analyze the effects of different recharge distributions on groundwater simulations. This study investigates the importance of the distribution of aquifer recharge in simulating regional groundwater flow in basin‐fill aquifers by calibrating a groundwater‐flow model to four different recharge distributions, all with the same total amount of recharge. Similarities are seen in results from steady‐state models for optimized hydraulic conductivity values, fit of simulated to observed hydraulic heads, and composite scaled sensitivities of conductivity parameter zones. Transient simulations with hypothetical storage properties and pumping rates produce similar capture rates and storage change results, but differences are noted in the rate of drawdown at some well locations owing to the differences in optimized hydraulic conductivity. Depending on whether the purpose of the groundwater model is to simulate changes in groundwater levels or changes in storage and capture, the distribution of aquifer recharge may or may not be of primary importance.     

Seismic dispersion and attenuation in layered porous rocks with fractures of varying orientations

Wave‐induced fluid flow plays an important role in affecting the seismic dispersion and attenuation of fractured porous rocks. While numerous theoretical models have been proposed for the seismic dispersion and attenuation in fractured porous rocks, most of them neglect the wave‐induced fluid flow resulting from the background anisotropy (e.g. the interlayer fluid flow between different layers) that can be normal in real reservoirs. Here, according to the theories of poroelasticity, we present an approach to study the frequency‐dependent seismic properties of more realistic and complicated rocks, i.e. horizontally and periodically layered porous rock with horizontal and randomly orienting fractures, respectively, distributed in one of the two periodical layers. The approach accounts for the dual effects of the wave‐induced fluid flow between the fractures and the background pores and between different layers (the interlayer fluid flow). Because C₃₃ (i.e., the modulus of the normally incident P‐wave) is directly related to the P‐wave velocity widely measured in the seismic exploration, and its comprehensive dispersion and attenuation are found to be most significant, we study mainly the effects of fracture properties and the stiffness contrast between the different layers on the seismic dispersion and attenuation of C₃₃. The results show that the increasing stiffness contrast enhances the interlayer fluid flow of the layered porous rocks with both horizontal and randomly orienting fractures and weakens the wave‐induced fluid flow between the fractures and the background pores, especially for the layered porous rock with horizontal fractures. The modelling results also demonstrate that for the considered rock construction, the increasing fracture density reduces the interlayer fluid flow while improves the dispersion and attenuation in the fracture‐relevant frequency band. Increasing fracture aspect ratio is found to reduce the dispersion and attenuation in the fracture‐relevant frequency band only, especially for the layered porous rock with horizontal fractures.     

The effect of orthorhombic anisotropy and its implication for oil recovery and reservoir exploitation

Results of an experimental study of shear‐ and compressional‐wave propagation in an orthorhombically anisotropic medium are presented. The experiments were performed on a physical model consisting of two sets of fractures. The first set consisted of orientated rubber inclusions simulating weak material‐filled cracks. The second set consisted of a system of closely spaced parallel fractures simulated by thin plates of epoxy resin, superimposed orthogonally on the first set. Three cases of fracture orientations within the model were identified and studied. Case 1 is analogous to a jointed fracture reservoir with one vertical set of fluid‐filled cracks or fractures and one non‐filled horizontal set. This case is referred to as JFV. Case 2 is analogous to a double fracture reservoir with one horizontal set of fluid‐filled fractures or cracks and one non‐filled vertical set. This is referred to as DFH. Case 3 is analogous to a double fracture reservoir with two vertical sets of fractures or cracks, with only one fluid‐filled. Case 3 is referred to as DFV. A pulse transmission method was performed on all three modelled cases along the three principal axes. A directional variation in the compressional‐ and shear‐wave velocities, as well as distinct shear‐wave splitting, was observed. The elastic constants for each case were determined and differences between them were noted and compared with the controlled results of both layered (transverse isotropy, TI) and vertically fractured (azimuthally anisotropic models, VF) media. The differences in elastic moduli and velocities indicate the potential of recognizing the different fracture orientations and suggest an approach to designing a method of drilling to further enhance oil recovery and reservoir exploitation.     

Roughness scale dependence of the relationship between tracer longitudinal dispersion and Peclet number in variable‐aperture fractures

The relationship between the longitudinal dispersion (D_L) and Peclet number (Pe) is crucial for predicting and simulating tracer through the variable‐aperture fracture. In this study, the roughness of the self‐affine fracture wall was decomposed into primary roughness (relatively large‐scale waviness) and secondary roughness (relatively small‐scale waviness) by a multiscaled wavelet analysis technique. Based on the complete dispersion mechanism (diffusion, macrodispersion, and Taylor dispersion) in the variable‐aperture fracture, three relationships (second‐order, power‐law, and linear relationships) between the D_L and Pe were investigated at large and small scales, respectively. Our results showed that the primary roughness mostly controlled the Taylor dispersion mechanism, whereas the secondary roughness was a dominant factor for the macrodispersion mechanism. Increasing the Hurst exponent and removing the secondary roughness led to the decreasing range of Pe where macrodispersion mechanism dominated the solute transport. It was found that estimating the D_L from the power‐law relationship based on Taylor dispersion theory resulted in considerable errors, even in the range of Pe where the Taylor dispersion mechanism dominated. The exponent of the power‐law relationship increased as the secondary roughness was removed. Analysing the linear relationship between the D_L and Pe revealed that the longitudinal dispersivity α_L increased linearly. However, this linear increase became weak as the Taylor dispersion mechanism dominated. In the range of Pe where the macrodispersion mechanism dominated, increasing the Hurst exponent caused the increase of α_L and the secondary roughness played a significant role in enhancing the α_L. As the Taylor dispersion mechanism dominated, the α_L was insensitive to the influence of multiscale roughness in variable‐aperture fractures.     

Seismic characterization of vertical fractures described as general linear-slip interfaces

Fluid flow in many hydrocarbon reservoirs is controlled by aligned fractures which make the medium anisotropic on the scale of seismic wavelength. Applying the linear‐slip theory, we investigate seismic signatures of the effective medium produced by a single set of ‘general’ vertical fractures embedded in a purely isotropic host rock. The generality of our fracture model means the allowance for coupling between the normal (to the fracture plane) stress and the tangential jump in displacement (and vice versa). Despite its low (triclinic) symmetry, the medium is described by just nine independent effective parameters and possesses several distinct features which help to identify the physical model and estimate the fracture compliances and background velocities. For example, the polarization vector of the vertically propagating fast shear wave S₁ and the semi‐major axis of the S₁‐wave normal‐moveout (NMO) ellipse from a horizontal reflector always point in the direction of the fracture strike. Moreover, for the S₁‐wave both the vertical velocity and the NMO velocity along the fractures are equal to the shear‐wave velocity in the host rock. Analysis of seismic signatures in the limit of small fracture weaknesses allows us to select the input data needed for unambiguous fracture characterization. The fracture and background parameters can be estimated using the NMO ellipses from horizontal reflectors and vertical velocities of P‐waves and two split S‐waves, combined with a portion of the P‐wave slowness surface reconstructed from multi‐azimuth walkaway vertical seismic profiling (VSP) data. The stability of the parameter‐estimation procedure is verified by performing non‐linear inversion based on the exact equations.     

Two-phase flow properties of a horizontal fracture: The effect of aperture distribution

A systematic numerical method has been presented to investigate the constitutive relationships between two-phase flow properties of horizontal fractures and aperture distributions. Based on fractal geometry, single rough-walled fractures are generated numerically by modified successive random addition (SRA) method and then aperture distributions with truncated Gaussian distribution are formed by shear displacement between lower and upper surfaces. (The truncated Gaussian distribution is used to describe aperture evolution under different normal stresses.) According to the assumption of two-dimensional porous media and local parallel plate model, invasion percolation approach is employed to model the two-phase flow displacement (imbibition) in generated horizontal fractures, in which capillary forces are dominant over viscous and gravity forces. For truncated Gaussian distributions, constitutive relationships from numerical simulation are compared to closed-form relationships and a good agreement is obtained. The simulation results indicate strong phase interference with the sum of two phase relative permeability values being less than one in the intermediate saturations. It is found that fracture properties related to residual saturations depend on spatial correlation of aperture distributions. Based on the simulation results, we proposed an empirical relationship between the fracture residual-saturation-rated parameters and the corresponding aperture distributions.     

Numerical examination of the factors controlling DNAPL migration through a single fracture

The migration of five dense nonaqueous phase liquids (DNAPLs) through a single fracture in a clay aquitard was numerically simulated with the use of a compositional simulator. The effects of fracture aperture, fracture dip, matrix porosity, and matrix organic carbon content on the migration of chlorobenzene, 1,2-dichloroethylene, trichloroethylene, tetra-chloroethylene, and 1,2-dibromoethane were examined. Boundary conditions were chosen such that DNAPL entry into the system was allowed to vary according to the stresses applied. The aperture is the most important factor of those studied controlling the migration rate of DNAPL through a single fracture embedded in a clay matrix. Loss of mass to the matrix through diffusion does not significantly retard the migration rate of the DNAPL, particularly in larger aperture fractures (e.g., 50 microm). With time, the ratio of diffusive loss to the matrix to DNAPL flux into the fracture approaches an asymptotic value lower than unity. The implication is that matrix diffusion cannot arrest the migration of DNAPL in a single fracture. The complex relationships between density, viscosity, and solubility that, to some extent, govern the migration of DNAPL through these systems prevent accurate predictions without the use of numerical models. The contamination potential of the migrating DNAPL is significantly increased through the transfer of mass to the matrix. The occurrence of opposite concentration gradients within the matrix can cause dissolved phase contamination to exist in the system for more than 1000 years after the DNAPL has been completely removed from the fracture.     

Evaluating morphological estimates of the aerodynamic roughness of debris covered glacier ice

Aerodynamic roughness length (z₀), the height above the ground surface at which the extrapolated horizontal wind velocity profile drops to zero, is one of the most poorly parameterised elements of the glacier surface energy balance equation. Microtopographic methods for estimating z₀ have become prominent in the literature in recent years, but are rarely validated against independent measures and are yet to be comprehensively analysed for scale or data resolution dependency. Here, we present the results of a field investigation conducted on the debris covered Khumbu Glacier during the post‐monsoon season of 2015. We focus on two sites. The first is characterised by gravels and cobbles supported by a fine sandy matrix. The second comprises cobbles and boulders separated by voids. Vertical profiles of wind speed recorded by a tower comprising five cup anemometers and deployed over both sites enable us to derive measurements of aerodynamic roughness that reflect their observed surface characteristics (0.0184 m and 0.0243 m, respectively). At the second site, z₀ also varied through time following snowfall (0.0055 m) and during its subsequent melt (0.0129 m), showing the importance of fine resolution topography for near‐surface airflow. To compare the wind profile data with microtopographic methods, we conducted structure from motion multi‐view stereo (SfM‐MVS) surveys across each patch and calculated z₀ using three previously published approaches. The fully three‐dimensional cloud‐based approach is shown to be most stable across different scales and these z₀ values are most correct in relative order when compared with the wind tower data. Popular profile‐based methods perform less well providing highly variable values across different scales and when using data of differing resolution. These findings hold relevance for all studies using microtopographic methods to estimate aerodynamic roughness lengths, including those in non‐glacial settings. Copyright © 2017 John Wiley & Sons, Ltd.     

Site-Specific Design for Dual Phase Recovery and Stabilization of Pooled DNAPL

Volume reduction and lowering of capillary pressure within a large DNAPL pool are utilized as objectives in the design of a large-scale dual phase recovery system at a chemical manufacturing facility in the United States. By reducing DNAPL pool height through mass removal, capillary pressure is lowered, resulting in a reduced potential for future vertical and horizontal mobilization of the chlorinated solvent DNAPL pool. The DNAPL pool extends over an approximately 200 m by 275 m area in low permeability fill deposits overlying a clay aquitard. A three-dimensional multiphase flow model was employed to arrive at a final design incorporating nine horizontal drains (total length 664 m) and a pulsed pumping system. The numerical model was calibrated to the results of a 42-day field pilot-test involving the removal of approximately 25,000 L of DNAPL from a single, 55 m long horizontal drain. Numerical simulation revealed that gravity drainage, as opposed to hydraulic gradients in the water phase, is the dominant recovery mechanism at this site. This stems from the relatively high density and the viscosity of the DNAPL, and the relatively low permeability of the formation deposits. The use of pulsed pumping is shown to reduce the volume of contaminated ground water recovered from the 9-drain system, without significant reduction of the total volume of DNAPL recovered.     

Bioaugmentation in a Well‐Characterized Fractured Rock DNAPL Source Area

A field demonstration was performed at Edwards Air Force Base to assess bioaugmentation for treatment of a well‐characterized tetrachloroethene (PCE) dense nonaqueous phase liquid (DNAPL) source area in fractured rock. Groundwater recirculation was employed to deliver remedial amendments, including bacteria, to facilitate reductive dechlorination and enhance DNAPL dissolution. An active treatment period of 9 months was followed by a 10‐month posttreatment rebound evaluation. Dechlorination daughter products were observed in both the shallow and deep fracture zones following treatment. In the shallow fracture zone, the calculated DNAPL mass removed was approximately equal to the DNAPL mass estimated using partitioning tracer testing, and no rebound in chlorinated ethenes or ethene was observed during the posttreatment period. A maximum DNAPL dissolution enhancement factor of 5 was observed in the shallow fracture zone. In the deep fracture zone, only approximately 45% of the DNAPL mass—as estimated via partitioning tracer testing—was removed and rebound in the total molar chlorinated ethenes + ethene was observed. The difference in behavior between the shallow and deep fracture zones was attributed to DNAPL architecture and the fracture flow field.     

Dependency of flow and transport properties on aperture distributions and compression states

Fluid conductivity and elastic properties in fractures depend on the aperture geometry – in particular, the roughness of fracture surfaces. In this study, we have characterized the surface roughness with a log-normal distribution and investigated the transport and flow behaviour of the fractures with varying roughness characteristics. Numerical flow and transport simulations have been performed on a single two-dimensional fracture surface, whose aperture geometry changes with different variances and correlation lengths in each realization. We have found that conventional measurement of hydraulic conductivity alone is insufficient to determine these two parameters. Transient transport measurements, such as the particle breakthrough time, provide additional constraints to the aperture distribution. Nonetheless, a unique solution to the fracture aperture distribution is still under-determined with both hydraulic conductivity and transport measurements. From numerical simulations at different compression states, we have found that the flow and transport measurements exhibit different rates of changes with respect to changes in compression. Therefore, the fracture aperture distribution could be further constrained by considering the flow and transport properties under various compression states.     

裂缝性储层的电各向异性响应特征研究

以水平裂缝分布的孔隙介质模型为基础,建立了水平裂缝面的平行裂缝分布和带粗糙表面的点接触裂缝模型,导出了两类裂缝模型的水平和垂直电阻率响应关系,分析了两类裂缝模型在没有围压和存在围压条件下的水平电阻率、垂直电阻率变化规律和电各向异性系数的变化特征.讨论了裂缝开度、裂缝密度和裂缝粗糙度等裂缝特征参数对裂缝性储层电各向异性的影响.为简化讨论,所有电性响应特征的分析都忽略了裂缝和孔隙表面的导电性和极化的影响.利用平行分布水平裂缝面模型和带粗糙面的裂缝模型,考察了围压条件下的裂缝性岩石的电各向异性响应特征,得到了对实际应用有意义的结果.     

Characterization of dipping fractures in a transverselyisotropic background

Although it is believed that natural fracture sets predominantly have near‐vertical orientation, oblique stresses and some other mechanisms may tilt fractures away from the vertical. Here, we examine an effective medium produced by a single system of obliquely dipping rotationally invariant fractures embedded in a transversely isotropic with a vertical symmetry axis (VTI) background rock. This model is monoclinic with a vertical symmetry plane that coincides with the dip plane of the fractures. Multicomponent seismic data acquired over such a medium possess several distinct features that make it possible to estimate the fracture orientation. For example, the vertically propagating fast shear wave (and the fast converted PS‐wave) is typically polarized in the direction of the fracture strike. The normal‐moveout (NMO) ellipses of horizontal reflection events are co‐orientated with the dip and strike directions of the fractures, which provides an independent estimate of the fracture azimuth. However, the polarization vector of the slow shear wave at vertical incidence does not lie in the horizontal plane – an unusual phenomenon that can be used to evaluate fracture dip. Also, for oblique fractures the shear‐wave splitting coefficient at vertical incidence becomes dependent on fracture infill (saturation). A complete medium‐characterization procedure includes estimating the fracture compliances and orientation (dip and azimuth), as well as the Thomsen parameters of the VTI background. We demonstrate that both the fracture and background parameters can be obtained from multicomponent wide‐azimuth data using the vertical velocities and NMO ellipses of PP‐waves and two split SS‐waves (or the traveltimes of PS‐waves) reflected from horizontal interfaces. Numerical tests corroborate the accuracy and stability of the inversion algorithm based on the exact expressions for the vertical and NMO velocities.     

Effects of fracture intersections on seismic dispersion: theoretical predictions versus numerical simulations

The detection and characterisation of domains of intersecting fractures are important goals in several disciplines of current interest, including exploration and production of unconventional reservoirs, nuclear waste storage, CO₂ sequestration, and groundwater hydrology, among others. The objective of this study is to propose a theoretical framework for quantifying the effects of fracture intersections on the frequency‐dependent elastic properties of fluid‐saturated porous and fractured rocks. Three characteristic frequency regimes for fluid pressure communication are identified. In the low‐frequency limit, fractures are in full pressure communication with the embedding porous matrix and with other fractures. Conversely, in the high‐frequency limit, fractures are hydraulically isolated from the matrix and from other fractures. At intermediate frequencies, fractures are hydraulically isolated from the matrix porosity but can be in hydraulic communication with each other, depending on whether fracture sets are intersecting. For each frequency regime, the effective stiffness coefficients are derived using the linear‐slip theory and anisotropic Gassmann equations. Explicit mathematical expressions for the two characteristic frequencies that separate the three frequency regimes are also determined. Theoretical predictions are then applied to two synthetic 2D samples, each containing two orthogonal fracture sets: one with and another without intersections. The resulting stiffness coefficients, Thomsen‐style anisotropy parameters, and the transition frequencies show good agreement with corresponding numerical simulations. The theoretical results are applicable not only to 2D but also to 3D fracture systems and are amenable to being employed in inversion schemes designed to characterise fracture systems.     

Controls on the aerodynamic roughness length and the grain‐size dependence of aeolian sediment transport

Estimates of the wind shear stress exerted on Earth's surface using the fully rough form of the law‐of‐the‐wall are a function of the aerodynamic roughness length, z₀. Accurate prediction of aeolian sediment transport rates, therefore, often requires accurate estimates of z₀. The value of z₀ is determined by the surface roughness and the saltation intensity, both of which can be highly dynamic. Here we report field measurements of z₀ values derived from velocity profiles measured over an evolving topography (i.e. sand ripples). The topography was measured by terrestrial laser scanning and the saltation intensity was measured using a disdrometer. By measuring the topographic evolution and saltation intensity simultaneously and using available formulae to estimate the topographic contribution to z₀, we isolated the contribution of saltation intensity to z₀ and document that this component dominates over the topographic component for all but the lowest shear velocities. Our measurements indicate that the increase in z₀ during periods of saltation is approximately one to two orders of magnitude greater than the increase attributed to microtopography (i.e. evolving sand ripples). Our results also reveal differences in transport as a function of grain size. Each grain‐size fraction exhibited a different dependence on shear velocity, with the saltation intensity of fine particles (diameters ranging from 0.125 to 0.25 mm) saturating and eventually decreasing at high shear velocities, which we interpret to be the result of a limitation in the supply of fine particles from the bed at high shear velocities due to bed armoring. Our findings improve knowledge of the controls on the aerodynamic roughness length and the grain‐size dependence of aeolian sediment transport. The results should contribute to the development of improved sediment transport and dust emission models. © 2018 John Wiley & Sons, Ltd.     

Analogue model studies of induction effects at auroral latitudes

In addition to field observations and numerical models, geomagnetic induction effects can be studied by scaled analogue model experiments. We present here results of analogue model studies of the auroral electrojet with an Earth model simulating the Arctic Ocean and inland conductivity structures in northern Fennoscandia. The main elements of the analogue model used were salt water simulating the host rock, an aluminium plate corresponding to the ocean and graphite pieces producing the inland highly conducting anomalies. The electrojet was a time-harmonic line current flowing at a (simulated) height of 100 km above northern Fennoscandia. The period simulated was 9 min.The analogue model results confirmed the well-known rapid increase of the vertical field when the coast is approached from the continent. The increase of the horizontal field due to induced ocean currents was demonstrated above the ocean, as well as the essentially negligible effect of these currents on the horizontal field on the continent.The behaviour of the magnetic field is explained with a simple two-dimensional thin-sheet model. The range, or the adjustment distance, of the ocean effect inland was found to be some hundreds of kilometers, which also agrees with earlier results of the Siebert-Kertz separation of IMAGE magnetometer data. The modelled inland anomalies evidently had too large conductivities, but on the other hand, their influence decayed on scales of only some tens of kilometers.Analogue model results, thin-sheet calculations, and field observations show that the induction effect on the horizontal magnetic field B_x near the electrojet is negligible. On the other hand, the vertical component B_z is clearly affected by induced currents in the ocean. Evidence of this is the shift of the zero point of B_z 0-1° southwards from the maximum of B_x. The importance of these results are discussed, emphasizing the determination of ionospheric currents.     

Numerical investigation of alternative fracture stiffness measures and their respective scaling behaviours

We study the mechanical deformation of fractures under normal stress, via tangent and specific fracture stiffnesses, for different length scales using numerical simulations and analytical insights. First, we revisit an equivalent elastic layer model that leads to two expressions: the tangent stiffness is the sum of an “intrinsic” stiffness and the normal stress, and the specific stiffness is the tangent stiffness divided by the fracture aperture at current stress. Second, we simulate the deformation of rough fractures using a boundary element method where fracture surfaces represented by elastic asperities on an elastic half‐space follow a self‐affine distribution. A large number of statistically identical “parent” fractures are generated, from which sub‐fractures of smaller dimensions are extracted. The self‐affine distribution implies that the stress‐free fracture aperture increases with fracture length with a power law in agreement with the chosen Hurst exponent. All simulated fractures exhibit an increase in the specific stiffness with stress and an average decrease with increase in length consistent with field observations. The simulated specific and tangent stiffnesses are well described by the equivalent layer model provided the “intrinsic” stiffness slightly decreases with fracture length following a power law. By combining numerical simulations and the analytical model, the effect of scale and stress on fracture stiffness measures can be easily separated using the concept of “intrinsic” stiffness. We learn that the primary reason for the variability in specific stiffness with length comes from the fact that the typical aperture of the self‐affine fractures itself scales with the length of the fractures.     

Numerical modelling of aeolian erosion over rough surfaces

The presence of non‐erodible roughness elements on erodible surfaces has the effect of absorbing part of the wind shear stress and thus protecting the erodible surface from wind erosion. This paper examines the shear stress distribution over roughness arrays of varying density, representing the progress of erosion on a bed of erodible and non‐erodible particles. Three‐dimensional numerical simulations, simulating wind flow over a bed of particles covered by roughness elements, were conducted in order to investigate the effect of roughness elements on the shear stress near the surface. The results of these simulations confirm that the erosion of soil by wind is strongly attenuated by the presence of roughness elements on the surface and depends on the geometric properties of the roughness elements. Based on the new numerical results obtained, a refinement of existing theoretical approaches is developed to describe the dependence of the friction velocity upon roughness frontal area and real exposed cover rate. The new formulation proposed will allow a more accurate evaluation of shear stress partitioning as a function of topographic changes. Copyright © 2010 John Wiley & Sons, Ltd.