Transmissivity Averaging in Fracture Flow on Self-affine Linear Profiles: Arithmetic, Harmonic, and Beyond

Numerical modeling of unidirectional flow in self-affine fractures using the lubrication approximation requires averaging of the transmissivity between the nodes. Seven averaging techniques are reviewed: arithmetic averaging of transmissivity; harmonic averaging of transmissivity; two averaging techniques derived by cell-based collocation method; global reconstruction of profile by means of multiquadrics; arithmetic averaging of aperture; harmonic averaging of aperture. In order to evaluate the performance of the seven techniques in terms of pressure errors and hydraulic aperture errors, self-affine profiles of 1024 nodes with different Hurst exponents (0.4 to 0.8) are generated. Every second node is then removed, resulting in 512-node profiles. Apertures at removed nodes are used in reference flow simulations on the 512-node profiles. Then, simulations with the seven averaging techniques are performed on 512-node profiles. Errors are computed with regard to the results obtained in the reference simulations. Reconstruction with multiquadrics is found to provide superior accuracy on self-affine profiles, followed by harmonic averaging of transmissivity or harmonic averaging of the aperture (Some of the errors analyzed in this study are minimized with the two last mentioned schemes.). Multiquadrics reconstruction is found to provide the best accuracy also on a smooth periodic profile.


Introduction
Fluid flow in the subsurface is often controlled by natural fractures.In order to predict flow and contaminant transport in such systems, hydraulic properties of individual fractures must be available.Rock fractures usually have rough walls which leads to a variable local fracture aperture, w, defined as distance between the walls.There is strong evidence that aperture distributions of fractures in different materials, including rocks, are self-affine, with the fractal dimension in the range of 2.2-2.8(Bouchaud et al. 1990;Måløy et al. 1992;Mandelbrot et al. 1984;Plouraboué et al. 1995;Schmittbuhl and Måløy 1997).Hydraulic properties of such fractures are described by hydraulic aperture.This is the aperture of a smooth-walled conduit that yields the same flow rate as the real, rough-walled fracture does, given the same pressure difference.Substantial effort has been invested over the past 3-4 decades into the study of hydraulic aperture of rough-walled fractures.A major part of these studies is numerical whereby a constant pressure difference is applied between two edges of a rectangular fracture or between two ends of a linear, one-dimensional fracture, and the resulting flow rate is computed by finite-difference or finite-element method.In the case of linear, one-dimensional fractures and under the assumptions of the lubrication theory approximation (Zimmerman et al. 1991), the numerical problem reduces to solving where L is the length of the fracture; P is the fluid pressure; c = c(x) is the local fracture transmissivity that can be defined as w 3 .Second-order finite-difference discretization of Eq. ( 1) on a uniform grid typically used in these models yields: where i is the node number.Both pressure and transmissivity are specified at nodes, which necessitates approximation, or averaging, of transmissivity, c, between the nodes, i.e., at locations i ± 1∕ 2 .Two types of transmissivity averaging are often used in numerical mod- els of subsurface flow: arithmetic and harmonic, given, respectively, by and Harmonic averaging is the method of choice in some hydrogeological simulators, e.g., MODFLOW (Langevin et al. 2017).It is often used in fracture flow models (Tsang and Neretnieks 1998) and porous media flow models (Aavatsmark 2002;Bessone et al. 2022;Roth et al. 1996;Colecchio et al. 2020).Its use for diffusion equation in strongly heterogeneous media was advocated by van Es et al. (2014).Its use in nonlinear diffusion problems was disputed by Kadioglu et al. (2008) who showed that it yields large errors when the transmissivity is P 3 .In fracture flow, harmonic averaging typically results in lower values of the hydraulic aperture than the arithmetic averaging because arithmetic average of two real numbers is always greater than or equal to their harmonic average.It was obtained by Lavrov (2022) that, in two-dimensional fracture flows, numerical schemes with arithmetic averaging perform better in terms of mass balance error.
Instead of averaging the transmissivity, some schemes make use of averaging the fracture aperture.Here, several options are available, too.The two main choices are arithmetic averaging and harmonic averaging.Arithmetic averaging of the local aperture was used, e.g., by Lenci et al. (2022a).It was pointed out by Lenci et al. (2022b) that harmonic averaging of aperture preserves the energy in the scheme. (1) The importance of averaging techniques in flux evaluation is not limited to fluid flow computations.Accurate evaluation of fluxes is crucial also for numerical heat transfer (Klepikova et al. 2021), solute transport (Moreno et al. 1988) and particulates transport (Lin et al. 2022;Mao et al. 2023) in fractures.
The objective of this work was to assess the performance of several averaging techniques in fracture flow on self-affine linear profiles.

Averaging Techniques
Seven averaging techniques were used in this study and are summarized in this Section.

Arithmetic Averaging of Transmissivity
Heuristically, arithmetic averaging of transmissivity is obtained by linear interpolation of transmissivity between the nodes, which yields Eq. ( 3).Another heuristic approach ending up with arithmetic averaging is by assuming the transmissivity to be constant between the nodes (Kadioglu et al. 2008); this implies discontinuous transmissivity at each node.
One way to properly derive Eq. ( 3) is by expanding the left-hand side of Eq. ( 1) as the derivative of a product: The derivation is provided in the Appendix of (Lavrov 2022) and is repeated here for reader's convenience.The aim is to derive a second-order accurate approximation of Eq. ( 5).A second-order approximation of the first term on the left-hand side of Eq. ( 5) by central difference is given by A second-order accurate approximation of the second term on the left-hand side of Eq. ( 5) can be obtained by adding Taylor expansions of c(x + Δx)P(x + Δx) and c(x-Δx)P(x-Δx), resulting in Using central difference approximations for d 2 P dx 2 and d 2 c dx 2 in Eq. ( 7) yields the following second-order approximation for the second term on the left-hand side of Eq. ( 5): Substitution of Eqs. ( 6) and (8) into Eq.( 5) recovers Eq. ( 2) with the internodal transmissivities given by Eq. ( 3): Thus, arithmetic averaging of transmissivity yields second-order accuracy as previously shown by Kadioglu et al. (2008) using a different technique.

Harmonic Averaging of Transmissivity
This type of averaging is hard to justify except heuristically.The latter can be done by imagining that the real transmissivity profile is replaced with a collection of conduits.Each conduit is centered around a grid node and has a constant aperture along its length (Fig. 1).
The fracture in Fig. 1 is assumed to be composed of three segments of constant aperture.The transmissivity of each segment is equal to the transmissivity specified in the node located in the middle of the segment, i.e., c i-1 , c i , c i+1 .Pressures are referred to the nodes as well, i.e., P i-1 , P i , P i+1 .Consider part of the fracture located within the shadow box.It consists of two half-segments connected in series.Pressure difference across these two half-segments is given by P i+1 -P i .Due to mass conservation, the flow rate through this part of the fracture is given by 2c i c i+1 c i +c i+1 P i −P i+1 Δx .Thus, c i±1∕ 2 is given by Eq. ( 4), with the following scheme as a result: By construction, this is a heuristic averaging method.The fact that it is based on the assumption of constant transmissivity within the cell has been recognized for at least two decades (Hyman et al. 2002).It can be shown that harmonic averaging yields second-order accuracy (Kadioglu et al. 2008), just as arithmetic averaging does.However, the magnitudes of the errors can be considerably different on some types of problems unless very fine grids are used (Kadioglu et al. 2008).
It should be noted that the assumption about piecewise-constant aperture profile implied in Eq. ( 4) is not realistic for fractures.Fractures have self-affine aperture maps.Thus, aperture varies at all scales, and it is only by chance that there could be a segment of constant aperture in a real fracture.
Both arithmetic and harmonic averaging procedures result in symmetric coefficient matrices.A symmetric matrix is necessary in order to ensure theoretical mass balance with the scheme (Romeu and Noetinger 1995).

Cell-Collocation-Based Averaging Techniques
When deriving the arithmetic averaging in Sect.2.1, another, alternative path could be taken when approximating the second term on the left-hand side of Eq. ( 5).Namely, we could have used central finite difference for each of the derivatives, dc/dx and dP/dx, resulting in the following scheme: (10) Fig. 1 Heuristic approach to harmonic averaging of transmissivity.The assumed shape of the fracture is shown with solid line.Nodes (i-1, i, i + 1) are located in the middle of constant-aperture segments i i+1 i-1

Δx
By construction, this is a second-order scheme.Another way to arrive at Eq. ( 11) is by using the cell-collocation method introduced in (Lau and Brebbia 1978).The use of cell-collocation method for deriving finite-difference approximations is detailed by Brebbia et al. (1984).Following the procedure outlined by Brebbia et al. (1984), we approximate P around node i by where the three basis functions are given by with the origin of the local coordinate system, x , being at node i.Since we are working with a three-node collocation cell, let us use the following interpolant for transmissivity between nodes i-1 and i + 1: The interpolant takes on the nodal values of c at the nodes.Substituting Eqs. ( 12) and (14) into Eq.(1), and using x = 0 as the collocation point (i.e., the point where the residual is equal to zero) yields: i.e., Inspection of Eq. ( 16) reveals that it implies the following averaging for the internodal transmissivities: We will call this collocation-based averaging.It is evident from Eq. ( 11) that the resulting coefficient matrix is not symmetric.As a consequence, this averaging technique yielded flux errors on the order of 1% in our simulations.Therefore, this averaging technique was not further considered in this study.
There are, however, at least two ways how to symmetrize Eq. ( 16) and thereby improve the flux error considerably.One way to symmetrize Eq. ( 16) is to use Eq. ( 17) with the lower sign for all internodal transmissivities except the one near the right boundary (see below for (11) the boundary treatment).Essentially, this means averaging c by means of a three-node stencil where one node is located to the left of the internodal point of interest, and two nodes are located to the right of the internodal point of interest.The resulting scheme follows immediately from Eq. ( 16) and is given by (The coefficient in front of P i ensures conservation.)We will call Eq. ( 18) right-collocation-based averaging in the remainder of this paper.The resulting coefficient matrix is symmetric.
Equation ( 18) has a four-node stencil (− 1, 0, 1, 2) and thus, cannot be used in the rightmost node.There, Eq. ( 16) should be used instead.This adjustment in one node retains the symmetry of the coefficient matrix.
Another way to symmetrize Eq. ( 16) is by using Eq. ( 17) with the upper sign for all internodal transmissivities except that near the left boundary.This means averaging c by means of a three-node stencil where two nodes are located to the left of the internodal point of interest, and one node is located to the right of the internodal point of interest.The resulting scheme is given by We will call Eq. ( 19) left-collocation-based averaging in the remainder of this paper.The resulting coefficient matrix is symmetric.
Equation ( 19) has a four-node stencil (− 2, − 1, 0, 1) and thus, cannot be used in the leftmost node.There, Eq. ( 16) should be used instead.This adjustment in one node retains the symmetry of the coefficient matrix.

Arithmetic Averaging of Aperture
Averaging methods considered so far in this Section were either heuristic (harmonic), or based on some formal mathematical derivation (arithmetic, collocation-based), or a combination of the two (right-and left-collocation-based).It should be remembered, however, that in the case of fracture flow the transmissivity is given by c = w 3 where w is the local aperture.Thus, another way to obtain the internodal transmissivities, c i±1/2 , is by reconstructing the fracture aperture profile between the nodes, which would yield the estimates of the local aperture, w i±1/2 .Such a reconstruction may be performed by either local or global approximation of the aperture profile.There is no unique way to construct such an approximation.In this study, three approximations are considered: two local (arithmetic and harmonic averaging of aperture) and a global approximation by multiquadrics.
Arithmetic averaging of aperture on a uniform grid is given by ( 18)

Harmonic Averaging of Aperture
Harmonic averaging of aperture is another way to perform local reconstruction.It is given by:

Global Reconstruction of Aperture Using Multiquadrics
Multiquadrics are a subclass of radial basis functions given, in general, by (Kansa 1990) where β ≥ -− /2; ε is the shape factor.In this study, multiquadrics with β = 1/2, originally introduced by Hardy (1971), were used.By using multiquadrics as interpolants, with the grid points serving as collocation points, the approximated aperture at i + 1/2 is given by (uniform grid with grid step Δx is assumed throughout this paper): where N is the number of grid points (N = 512 in this study);  = εΔx is the dimension- less shape factor.The condition number of the linear system used to obtain the interpolation coefficients, C k , increases rapidly as the shape factor, ε, decreases.As ε increases, multiquadrics become more 'pointy', and the interpolation Eq. ( 23) approaches the arithmetic averaging of aperture, Eq. ( 20).Various strategies are available for choosing ε.For instance, different (e.g., random) individual shape factors can be assigned to different k (Sarra,Sturgill 2009).In this study, the simplest possible choice was made by using the same shape factor, ε = 1.0, for all quadrics.With this choice of ε, and N = 512, the condition number is sufficiently low (2.1•10 6 ).

Methodology
None of the averaging techniques presented in Sect. 2 appear to be preferable based on the derivation only.It should be remembered though that Eq. ( 1) is to be solved on selfaffine linear profiles.In order to evaluate which averaging procedure suits best for solving Eq. ( 1) on such profiles, the following approach was employed in this study: First, a profile of 1024 nodes was constructed.Then, every second node was removed from the 1024-node profile.The hydraulic aperture of the resulting 512-node profile was computed, using the seven types of averaging introduced in Sect. 2 (arithmetic transmissivity, harmonic transmissivity, left-and right-collocation-based, reconstruction by global multiquadric approximation, arithmetic aperture, and harmonic aperture).The results of the seven computations were compared to the 'true' (reference) value which was obtained with the internodal transmissivities corresponding to the 'true' internodal apertures, i.e., the apertures at the removed nodes of the 1024-node profile.Profiles of 1024 nodes each were constructed as follows: Five square self-affine aperture maps of 1024 × 1024 nodes were generated using the spectral synthesis method (Huang et al. 1992).The Hurst exponent of these surfaces was H = 0.4, 0.5, 0.6, 0.7 and 0.8.Lower Hurst exponent means that the map appears rougher.It is currently believed that most fractures have H = 0.8 which is sometimes considered a "universal value" (Lenci et al. 2022b).However, lower values of H, down to 0.45, have been measured as well (Lenci et al. 2022b).The value of the Hurst exponent set in the spectral generator was verified for each surface using a recently introduced version of box-counting method (Ai et al. 2014) and was double-checked with the variogram method (Ruello et al. 2006).From each of these surfaces, 10 linear profiles of 1024 nodes were extracted in the x-direction and 10 linear profiles of 1024 nodes were extracted in the y-direction.Each of these linear profiles was offset to make its average height equal to zero.The resulting profile was rescaled so that its range be equal to a given value, Δw. ('Range' here refers to the difference between the highest peak and the deepest valley.)The following values of Δw were used: 1, 2, 3, 4 and 5 mm.In the rescaled profile, the minimum height was identified.The rescaled profile was then offset so that the minimum height be equal to a given value, w min , where w min was 1, 2, 3, 4 and 5 mm.The ratio of root mean square variation (RMS) to the mean aperture varied in the range 0.0278 to 0.3091 for x-profiles and 0.0257 to 0.3008 for y-profiles.In each of these rescaled-and-offset profiles, every other node was removed thus giving rise to a 512-node profile.The pressure diffusion equation was solved on this profile with the seven averaging techniques discussed above.The pressure applied at the left-hand boundary was equal to 1, the pressure applied at the right-hand boundary was 0. In addition, on each 512-node profile, a reference simulation was performed with the 'true' internodal transmissivities based on the apertures at the removed nodes.The total number of flow computations excluding the reference simulations was thus equal to 17,500: 5 Hurst exponents × 20 linear profiles (10 in x-and 10 in y-direction) × 5 range values × 5 offset values × 7 averaging schemes.Direct solvers were used in flow computations and in solving for multiquadric coefficients.
In each simulation, the mass balance error was evaluated as follows: The flux between the two leftmost nodes, Q l , and the two rightmost nodes, Q r , was evaluated.The mass balance error was the magnitude of the difference between the two divided by their average, It was required for the mass balance error to be below 10 -6 for a simu- lation to be accepted.
In each simulation, the error in hydraulic aperture and the error's absolute magnitude were calculated.The error is given by w (a)  h − w (t) h w (t)  h where a stands for averaging type (arithmetic transmissivity, harmonic transmissivity, left-collocation, right-collocation, multiquadrics, arithmetic aperture, harmonic aperture); w (t)  h is the 'true' value of the hydraulic aperture obtained in the reference simulation.
For each value of the Hurst exponent and each type of averaging scheme, average error in hydraulic aperture and average absolute value of the error in hydraulic aperture were calculated by summing the individual values obtained in each simulation and dividing by the number of simulations (250, i.e., 10 profiles × 5 range values × 5 offset values).This was done separately for profiles extracted in the x-and y-directions.
In each simulation, three norms were calculated to assess the error in the pressure distribution along the profile: where summation was performed over internal nodes of the grid (510 in total); the meaning of superscripts is explained above. 1or each value of the Hurst exponent and each type of averaging scheme, average norms were calculated by summing the individual values obtained in each simulation and dividing by the number of simulations (250).This was done separately for profiles extracted in the x-and y-directions.
In the simulations described so far in this Section, the minimum aperture, w min , was 1, 2, 3, 4 or 5 mm and the aperture range, Δw, was 1, 2, 3, 4 or 5 mm.In order to investigate the performance of the seven schemes with smaller apertures, an additional series of simulations was conducted, with w min = 0.1 mm and Δw = 1, 2, 3, 4 or 5 mm.Only the 1024 × 1024 fracture with Hurst exponent equal to 0.8 was used to construct linear profiles as this H value is believed to be the most common for rock fractures (Lenci et al. 2022b).As before, 10 profiles were selected in x-direction and 10 in y-direction.

Results Obtained with w min = 1 mm to 5 mm
The average errors obtained with profiles extracted in the x-and y-directions are displayed in Figs. 2 and 3, respectively.It follows from these Figures that, in the case of x-profiles, the best results are obtained with multiquadric approximation of the aperture or with harmonic averaging of the aperture.The errors obtained with these two techniques are virtually indistinguishable from each other, no matter what type of error is used as a selection criterion.
Average errors obtained with left-and right-collocation-based averaging are close, as evidenced in Figs. 2 and 3.This is hardly surprising since the schemes are similar.
Examining the errors obtained in individual simulations revealed that the magnitude of error in hydraulic aperture obtained with arithmetic averaging of transmissivity is higher than that obtained with harmonic averaging of transmissivity in the majority of cases, especially for lower H (Table 1).This is also evident from Figs. 2 and 3. Thus, if the choice is to be made between arithmetic averaging of transmissivity and harmonic averaging of transmissivity, the latter appears for self-affine linear profiles to be a better choice. (24) Errors obtained with arithmetic averaging of transmissivity and displayed in Figs. 2 and 3 are relatively large.In order to understand the poor performance of arithmetic averaging, it is instructive to take a closer look at individual errors in the simulations rather than at their averaged values.In Fig. 4, errors in hydraulic aperture obtained with arithmetic averaging of transmissivity and with multiquadric reconstruction (the best performing scheme) are displayed as a function of 'true' hydraulic aperture for 250 simulations obtained on H = 0.4 x-profiles.Each point represents an individual simulation.As evident from Fig. 4, errors obtained with multiquadric reconstruction are distributed symmetrically around 0, while errors obtained with arithmetic averaging of transmissivity are predominantly positive.It means that arithmetic averaging of transmissivity systematically overestimates the hydraulic aperture.A graph similar to the one in Fig. 4 is obtained also if, instead of multiquadrics reconstruction, harmonic average of transmissivity is used.

Results Obtained with w min = 0.1 mm
All average errors, except signed error in the hydraulic aperture, were at their minimum when multiquadric reconstruction was used.In particular, all pressure errors (L 1 , L 2 , L ∞ ) were minimized this way.Only the average signed error in hydraulic aperture was minimized with an averaging technique other than multiquadrics, namely with harmonic averaging of transmissivity.For all other types of errors, the latter was the second best averaging technique.This is in contrast to the simulations discussed in Sect.4.1 where the second-best technique was typically harmonic averaging of the aperture.
Thus, three averaging techniques stand out for flow simulations on 1D self-affine fracture profiles: (1) multiquadrics reconstruction, (2) harmonic averaging of transmissivity and (3) harmonic averaging of aperture.Multiquadrics reconstruction of the aperture profile appears, on average, to be superior, followed by the two others.The higher accuracy of the multiquadrics reconstruction is bought by the additional computational cost this technique entails, compared to all other techniques.

Discussion
Multiquadrics have, in the past, demonstrated surprisingly good performance in various interpolation and reconstruction jobs even though it was pointed out that the reasons for this good performance were unknown (Kansa 1990).Our work provides yet another example of an application where reconstruction with multiquadrics can be used with advantage, compared to other, traditionally used, techniques, e.g., arithmetic or harmonic averaging of transmissivity.
The choice of the shape factor for multiquadrics in this work was based solely on the requirement of sufficiently small condition number.No attempt was made to optimize further.
Harmonic averaging of aperture or transmissivity yields virtually the same accuracy as the multiquadrics reconstruction.Computational overhead of the latter is, however, higher.The fact that the second best option is provided by harmonic averaging schemes is consistent with the earlier results by Romeu and Noetinger (1995) who found harmonic averaging of transmissivity to be a superior scheme for linear flow in porous media.They did not consider multiquadrics reconstruction though (the multiquadrics-based method makes sense for fracture flow where it has a clear geometric meaning).
The focus in this paper has been on the flow along random profiles.It is instructive to compare the performance of the seven averaging techniques on a smooth profile for which an analytical solution is available.As an example, consider transmissivity given by ( 27) Equation ( 27) represents a smooth periodic transmissivity.The analytical solution to Eq. ( 1) with P 0 = 1, L = 1 and c(x) given by Eq. ( 27) is The analytical solution is displayed in Fig. 5 along with the numerical solution obtained on a 512-node grid with multiquadric reconstruction of c 1/3 .Errors in hydraulic aperture obtained with the seven averaging schemes and N = 512 nodes are summarized in Table 2.
From Table 2, multiquadrics approximation yields, with a good margin, the best accuracy on a smooth periodic profile.This is contrary to the results obtained with self-affine profiles where multiquadrics reconstruction and harmonic averaging of the aperture produced very similar errors.This is yet another illustration that there is apparently no universal choice of the averaging technique and the choice should be determined by the nature of the problem at hand, a conclusion made by Kadioglu et al. (2008).
Excellent performance of multiquadrics on a profile given by Eq. ( 27) should perhaps come as no surprise since the profile is smooth and the global reconstruction technique does a superior job here (at the expense of some computational overhead).
It should be noted that we called multiquadrics approximation, arithmetic averaging of aperture and harmonic averaging of aperture "reconstruction techniques" in this text.However, all seven averaging techniques examined in this study could be considered as reconstruction techniques since in all of them we attempt to assign transmissivity to a location between the nodes.
Uniform grids were used in this study.Uniform grids are typically used in fracture flow simulations.Averaging on nonuniform grids should be investigated in the future.
Only linear (1D) fracture profiles were considered in this study.In practice, fracture flow simulations are performed on two-dimensional aperture landscapes.Dimensionality has a crucial effect on conductivity and flow structure, as demonstrated, e.g., by Colecchio et al.(2021).The results of our study might therefore not be directly applicable to 2D grids.In particular, in 1D, even a single node with zero aperture results in a trivial solution (zero flow rate), and a simulation does not need to be performed.Treatment of zero apertures is therefore not relevant for 1D schemes.In 2D, zero apertures may still result in a conductive fracture, and a simulation can be performed.In this case, in order to avoid a singular coefficient matrix, zero apertures are usually replaced with very small values.Comparison of different averaging schemes on 2D grids, especially with clusters of nodes having very small apertures, should be the subject of a future work.It should be noted, however, that it might be difficult to match the number of 1D simulations used in this study (a total of 17,500) when doing 2D simulations.
The work presented in this paper focused on the effects that averaging has on the flow.It might be worthwhile, in the future, to examine also how these seven reconstruction techniques affect geometric properties of the reconstructed fractures, e.g., their correlation length, RMS, Hurst exponent, etc.
All development in this paper has been for Newtonian fluids.The situation would become more complicated with non-Newtonian fluids.This can be illustrated using the ( 28) simplest possible non-Newtonian rheology, i.e., a power-law fluid.For a linear fracture flow under the lubrication theory approximation, the fluid velocity averaged across the fracture aperture is given for such a fluid by equation (Lavrov 2023) where n is the flow index; C is the consistency index.Mass conservation for one-dimensional steady incompressible flow, i.e., ∇(w ) = 0 , results in a p-Laplace equation: Inspection of Eqs.(29-30) suggests that the very concept of local fracture conductivity as a proportionality factor between velocity and pressure gradient acquires an additional dimension here: it becomes dependent on the pressure gradient.One promising approach here could be using reconstruction of aperture instead of averaging of transmissivity.This approach was pursued, e.g., in the recent work by Lenci et al.where arithmetic averaging of w between nodes was used to approximate the internodal aperture (Lenci et al. 2022b).As another alternative, the method of aperture reconstruction with multiquadrics advocated in the present paper could be tried for Eq.(30).Whether this method has advantages as compared to other averaging techniques also in the case of non-Newtonian fluids should be investigated in a dedicated study.
It should be stressed that conclusions drawn in our study pertain to self-affine linear profiles.It is straightforward to construct an example of profile where local averaging techniques, e.g., arithmetic or harmonic averaging of aperture or transmissivity, should perform at least as well as multiquadrics.An example I will use here is inspired by the recent work on 'binary media' (Colecchio et al. 2021) where fractures were constructed such that local aperture could take only two values,: w l or w u .Let us construct a special version of such 'binary media', namely a linear profile of 100 nodes where all nodes with index i < 50 have w = w l , while all nodes with i > 49 have w = w u .On such a profile, any local averaging method, e.g., arithmetic or harmonic averaging of aperture or transmissivity, should perform at least as well as any global reconstruction technique such as multiquadrics reconstruction.
Moreover, the very methodology employed in our study, with removing every other node and using the original profile for a reference computation, was designed for selfaffine profiles and may become faulty in other cases, e.g., when the aperture profile has a distinctive regular structure.As an illustration, let us again construct a 'binary' profile of a 100 nodes where every node with an even index has w = w l , while every node with an odd index has w = w u .This resembles the profile studied, e.g., by Haugerud et al. (2022).This is the profile to be used for the reference simulation.Now, if every other node is removed from this profile and a flow simulation is performed, all seven averaging methods considered in this study will perform poorly.From the averaging point of view, there is no way to distinguish between two possible original profiles once every other node has been removed: the one where all apertures were the same and the one where even nodes had w = w l and odd nodes had w = w u .Once the nodes have been removed, the information is lost (cf.Nyquist theorem).This brings us to the next issue, that of correlation length in the original profile.In practice, various degree of correlation may exist in an aperture profile, with the correlation length being determined by various factors, e.g., by the magnitude of the shear displacement (Brown et al. 1986;Wang et al. 1988).For instance, if two identical (matching) self-affine fracture surfaces are put together with zero offset, the resulting aperture profile will be w = const (i.e., infinite correlation length), and any local averaging technique in a numerical flow simulation will perform perfectly [cf.also the results of Nicholl et al. (1999) where several local averaging techniques were tested and their results converged at high mesh resolution, namely when the grid spacing became smaller than 1/5 of the correlation length].Global reconstruction techniques, e.g., based on multiquadrics, would have no advantage in this case.Thus, there might be an extra angle to the problem of choosing between different averaging techniques in fracture flow modeling.The role of correlation length and profile structure in this choice should be examined in the future.Finally, it should be emphasized that this study makes use of the lubrication theory approximation when formulating the flow problem.Our findings suggest that multiquadrics  reconstruction provides superior results when the flow is modeled using the lubrication approximation.This does not suggest that multiquadrics reconstruction provides the most accurate flow solution in general.The effect of out-of-plane flow on the accuracy of 1D flow simulations has been recognized in the research community for at least 25 years by now (Oron and Berkowitz 1998;Nicholl and Detwiler 2001;Brush and Thomson 2003;Wang et al. 2015Wang et al. , 2018Wang et al. , 2020;;He et al. 2021).An earlier attempt to reconcile the use of Eq. ( 1) with experimental results obtained using Hele-Shaw cells with artificial rough surfaces (not self-affine), while solving the flow with different averaging techniques, can be found in the work by Nicholl et al. (1999) There, corrections were introduced in order to account for out-of-plane flow dynamics, all with limited success.

Conclusions
The following conclusions are drawn: • In most cases studied here, global reconstruction of the aperture profile using multiquadrics provided superior accuracy in terms of both the modeled hydraulic aperture and the pressure error for fluid flow on one-dimensional self-affine aperture profiles.• The second best averaging technique was either harmonic averaging of the aperture or harmonic averaging of transmissivity.The former can outperform multiquadrics with regard to certain types of errors when the ratio of the minimum aperture to the aperture range is sufficiently large.The latter can outperform multiquadrics reconstruction with regard to certain types of errors when the ratio of the minimum aperture to the aperture range is sufficiently small.• On a smooth periodic profile, multiquadrics reconstruction outperformed all other techniques examined here, with a good margin.

Fig. 2
Fig.2Errors averaged over 8750 simulations performed on 512-node profiles oriented in the x-direction: a absolute value of error in hydraulic aperture; b error in hydraulic aperture; c pressure solution error using L 1 ; d pressure solution error using L 2 ; e pressure solution error using L ∞ .Results obtained on 1024-node profiles were used as reference when computing errors in each 512-node simulation ▸ Fig. 2 (continued)

Fig. 3
Fig.3Errors averaged over 8750 simulations performed on 512-node profiles oriented in the y-direction: a absolute value of error in hydraulic aperture; b error in hydraulic aperture; c pressure solution error using L 1 ; d pressure solution error using L 2 ; e pressure solution error using L ∞ .Results obtained on 1024-node profiles were used as reference when computing errors in each 512-node simulation ▸ Fig. 3 (continued)

Fig. 4 Fig. 5
Fig. 4 Signed error in hydraulic aperture obtained in individual simulations with arithmetic averaging of transmissivity and with multiquadric reconstruction vs. 'true' hydraulic aperture.Data for simulations with H = 0.4 are displayed.Each point represents an individual simulation

Table 1
Fraction of simulations in which the magnitude of error in hydraulic aperture obtained with arithmetic averaging of transmissivity was greater than with harmonic averaging of transmissivity