Bayesian full-waveform inversion of tube waves to estimate fracture aperture and compliance

. The hydraulic and mechanical characterization of fractures is crucial for a wide range of pertinent applications, such as, for example, geothermal energy production, hydrocarbon exploration, CO 2 -sequestration, and nuclear waste disposal. Direct hydraulic and mechanical testing of individual fractures along boreholes does, however, tend to be slow and cumbersome. To alleviate this problem, we propose to estimate the effective hydraulic aperture and the mechanical compliance of isolated fractures intersecting a borehole through a Bayesian Markov chain Monte Carlo (MCMC) inversion of full-waveform tube-5 wave data recorded in a vertical seismic proﬁling (VSP) setting. The solution of the corresponding forward problem is based on a recently developed semi-analytical solution. This inversion approach has been tested for and veriﬁed on a wide range of synthetic scenarios. Here, we present the results of its application to observed hydrophone VSP data acquired along a borehole in the underground Grimsel Test Site in the Central Swiss Alps. While the results are consistent with the corresponding evidence from televiewer data and exemplarily illustrate the advantages of using a computationally expensive stochastic, instead of a 10 deterministic, inversion approach, they also reveal the inherent limitation of the underlying semi-analytical forward solver.

such as groundwater management, geothermal energy production, hydrocarbon exploration, CO 2 -sequestration, and nuclear waste disposal. If a tube wave is generated at a fracture due to an incident P-wave, the amplitude ratio of the two wave types can be used to estimate fracture compliances (e.g., Bakku et al., 2013) or fracture permeability (e.g., Hardin and Toksöz, 1985;Li et al., 1994), while the amplitude ratio of the P-wave-induced tube waves to the S-wave-induced tube waves can be inverted for the orientation of fractures (e.g., Lee and Toksöz, 1995). The algorithm of Hornby et al. (1989) uses the arrival times of 40 reflected tube waves to invert for the locations of permeable fractures and the reflectivity of tube waves to estimate the effective aperture of fractures. In the field of seismoelectrics, Zhu et al. (2008) showed that tube waves create electromagnetic waves when encountering fractures, which also have the potential to be used for fracture characterization.
The above methods do, however, require extensive manual conditioning of the data, like amplitude picking or time-gating of events. Furthermore, they are unable to provide an estimate of uncertainty and/or to identify multiple solutions that are 45 equally likely. The objective of this work is to alleviate these limitations by providing an algorithm that considers the entire wavefield for characterizing fractures in terms of their hydraulic apertures and mechanical compliances as well as the associated uncertainties with a minimal amount of human interaction. To this end, we propose a Bayesian full-waveform inversion approach in combination with a recent semi-analytical approach  as an efficient and robust forward solver. The proposed algorithm uses as input the complete P-and tube-wave fields with minimal prepro-50 cessing to invert for the effective hydraulic fracture aperture, the mechanical fracture compliance, the bulk-and shear-modulus of the background rock as well as some auxiliary parameters. We use a stochastic inversion algorithm in order (1) to obtain an entire ensemble of solutions, which, in turn, provides a measure of uncertainty and (2) to account for the strong non-linearity of the problem and to avoid getting stuck in local minima. We first present our stochastic full-waveform inversion approach, followed by a synthetic example and an application to field data from the underground Grimsel Test Site (www.grimsel.ch) in 55 Switzerland and a subsequent discussion of the results.

Method
The goal of our stochastic inversion approach is to estimate the posterior probability density function (PDF) p(m|d), which in stochastic terms describes the adequacy of a model m given the observed data d. We do this by relying on the following approximation of Bayes' theorem (Bayes, 1763): where p(m) is the prior PDF describing any a priori knowledge we have about the model parameters and L(m|d) is the likelihood quantifying how well a model m explains the data d. Following Tarantola (2005), we define the likelihood as: where D and σ e are the amount of data points and the standard deviation of the data-error, respectively. The forward operator G calculates synthetic data d syn based on a model m: We use a novel semi-analytical algorithm for G, which evaluates the Green's function analytically in the frequency-space domain for a zero-offset VSP setting . This is done in parallel for a limited number of individual frequencies. Then, Green's functions for the complete frequency band are obtained by spline interpolation. The frequencies,

70
for which the Green's functions are actually calculated, are selected such that the maximum error caused by the interpolation (i.e., the difference between an interpolated and a fully calculated dataset) is two orders-of-magnitude smaller than the largest value in the dataset. After multiplication with the Fourier transform of the source wavelet and a subsequent inverse Fourier transformation, we obtain the full-waveform signals in the time-space domain.
In the considered forward operator G, seismic tube waves are generated and scattered at fractures characterized by their static 75 apertures L 0 and compliances Z. A tube wave is generated when a P-wave hits a fluid-filled fracture intersecting the borehole, as the fracture is compressed and fluid is injected into the borehole. We describe this process in the frequency domain for a horizontal fracture with the tube-wave generation potential φ g : where N is the number of fractures in the medium, ρ f the density of the fluid and δ the Dirac delta function. Depth is denoted 80 by z and sub-or superscripts i refer to the ith fracture. Note that this formulation requires the depth vector z to explicitly sample the depth levels of all prevailing fractures. Therefore, the sampling along z determines the minimal distance between two adjacent fractures that can be resolved. The tube-wave velocity c T is given by (White, 1983): with K f and µ being the fluid bulk modulus and the shear modulus of the formation, respectively. The pressure fields of the 85 tube wave p (i) t and the incoming P-wave p (i) inc are then given by: where σ 0 is the amplitude of the normally incident plane P-wave, j = √ −1 the imaginary unit, ω the angular frequency, k r the radial wavenumber for a rigid, non-deformable fracture (a function of L 0 ), α f the fluid velocity, α eff the effective fluid velocity 90 in the fracture (a function of L 0 and Z), and R the borehole radius. H n denotes the Hankel function of the first kind of order n, ζ the effective radial wavenumber (a function of L 0 and Z) and ρ the density of the embedding background rock. V P and V S are the P-wave and S-wave velocity in the background rock, respectively. Note that σ 0 drops out of equation 4 due to the ratio of p inc . When a tube wave propagating along the borehole interface encounters a fracture, fluid flow from the borehole into the 95 fracture is triggered. This leads to reflection and transmission of tube waves. This process is described with the scattering potential φ s in the frequency domain: where η is the interface compliance given by:

100
Note that the interface compliance differs from the fracture compliance. It linearly relates the velocity discontinuity ∆V across the fracture to the acoustic pressure p: ∆V = jωηp . Note also, that in our implementation of this forward solver, tube waves that are generated at borehole enlargements, such as, for example, washouts and bit-size changes, or at high acoustic impedance contrasts due to lithological changes are not taken into account. Further details about the tube-wave generation and scattering potentials, and the algorithm itself, can be found in .

105
For the forward operator G as described so far, we assumed the fractures to be horizontally oriented. To account for arbitrary incidence angles we have extended the above algorithm for the forward operator G, following the description given by .
To improve the estimation of the fracture compliance Z we have extended the forward operator of  to include transmission losses of P-waves across fractures, by using the angle-dependent transmission coefficient described by 110 the linear slip theory (Schoenberg, 1980). Accordingly, the P-and S-wave reflection coefficients R P and R S , as well as the Pand S-wave transmission coefficients T P and T S , for an incoming P-wave are given by: where with the superscript m being 1 for the medium above and 2 for the medium below the fracture. The angles θ m and ψ m refer to the P-wave and the S-wave reflection angles if the superscript m is 1, and the corresponding transmission angles if the superscript m is 2. Z T , Z N , and ρ denote the fracture compliance in the transverse direction (parallel to the fracture), the 120 fracture compliance in the normal direction (perpendicular to the fracture) and the density, respectively. Note that in this study, we assume for simplicity that Z = Z T = Z N . We solve equation 10 for the four coefficients, but we only use the transmission coefficient T P to reduce the amplitude of the P-wave after having crossed a fracture, because we do not consider reflections or S-waves in this study.
In order to fit the observed data, we implemented the forward operator of  such that the following 125 features are explicitly included: (1) Geometrical spreading of P-waves is accounted for by multiplying equation 7 with 1/z.
Note that other attenuation mechanisms of the P-wave, besides geometrical spreading and transmission losses across fractures, are neglected.
(2) The algorithm assumes a uniform embedding background medium. To accommodate for P-wave velocity changes above the considered borehole section, we introduce a variable source depth. This is an auxiliary parameter estimated during the inversion. (3) The algorithm assumes an isotropic background medium. As the particle motion of a P-wave is 130 different compared to that of a tube wave in the elastic medium surrounding the borehole, the two wave types are sensitive to the background medium properties in different directions. Therefore, taking anisotropy into account is important for fitting observed data. We do this by estimating different effective isotropic shear moduli for the P-wave and for the tube wave. Thus, the shear modulus µ in equation 5 becomes µ t , the tube-wave shear modulus.
Due to the non-linearity of the problem, we cannot infer the posterior PDF directly, but need to infer it by sampling the prior 135 PDF and the likelihood according to relation 1. For this, we chose to use a Markov chain Monte Carlo (MCMC) approach. This algorithm walks randomly through the solution space accepting or rejecting proposed models m prop , which are drawn from a symmetric proposal distribution, with the Metropolis acceptance probability α (Metropolis et al., 1953): towards the posterior PDF due to parallel and interacting Markov chains as well as a model-proposal scheme that uses a database of previously accepted models to avoid sampling areas of low posterior probability and focusing on the interesting areas of the solution space.
The viability and accuracy of the algorithm have been tested and verified on a variety of synthetic case studies, an example 145 of which is shown in the next section. Subsequently, we apply the proposed inversion scheme to hydrophone VSP data acquired at the underground Grimsel Test Site in the Central Swiss Alps.
3 Results: A synthetic example with real noise Before applying our inversion algorithm to observed data, we have run tests on synthetic data to ensure that the algorithm performs as expected. As in these experiments the same forward solver was used for the generation and the inversion of the 150 data, the corresponding results only allow to draw conclusions with regard to the inversion algorithm itself, but not with regard to the information content of the data. The test case shown here features two fractures at 10 and 19 m depth. The receiver spacing is 1 m. To make this synthetic study more pertinent and challenging, we contaminated the dataset with actual ambient noise from a corresponding field dataset at the underground Grimsel Test Site in Switzerland. The resulting data are plotted in Figure 1a. This synthetic test differs from the field-data example shown in the next section in two ways: (1) It uses as a forward solver the algorithm proposed by  and  without taking transmission losses, geometrical spreading for P-waves, velocity changes above the considered borehole section, or anisotropy into account, because these features are not present in the underlying synthetic data.
(2) While the wavelet is based on a mean trace for the field data, we treat it as unknown and, thus, estimate it in the synthetic example. We do this, by inferring the coordinates of six pilot points, 160 from which we obtain the wavelet by a shape-preserving piecewise cubic interpolation (Hunziker et al., 2019).
The inversion was run once with three parallel Markov chains. Figure 2 shows the estimate of the hydraulic fracture aperture and the mechanical compliance for the two fractures as a function of the number of forward simulation steps. For all four unknowns, the three chains converge nicely to the true values. This behavior illustrates, that the algorithm works properly even when the data are contaminated with correlated, realistic noise. Simulated data based on a model proposed at the end of the first Markov chain agree very well with the input data ( Figures   1a, b). Note that besides the direct P-wave () and the tube waves generated at fractures (), the tube waves reflected at fractures () are also visible. The latter are visible neither in the noise-contaminated input data nor in the actual field data.
The VSP data, considered in the following, were recorded in crystalline rocks at the underground Grimsel Test Site in Switzer-170 land using a 12-receiver hydrophone string with a receiver spacing of 1 m. In the course of the experiment, the hydrophone string was repositioned, such that, the recorded traces are separated by 0.5 m. The borehole had a diameter of 0.147 m. As a source, a single-handed 2 kg hammer was used at the wellhead, which excited frequencies between 0.1 and 4 kHz with a peak around 1.5 kHz. In this study, we consider a 20-m-long subsection between 17 and 37 m depth, consisting of 41 hydrophone receiver positions. Through visual inspection of the VSP dataset, complemented by evidence from optical and acoustic televiewer 175 data (Krietsch et al., 2018), three fractures at 23.5, 23.9 and 25 m depth have been identified.
Preprocessing of the data included a gentle bandpass filter to suppress high-frequency noise, a static shift correction to remove positioning errors, and a cosine taper to blank out the later arriving S-wave and associated tube waves. The data after preprocessing are shown in Figure 3. The P-wave and the tube waves are clearly visible. However, scattered tube waves, as described by equation 8, are weak in amplitude and drop below the noise level. As the first and the second fracture are located 180 closely together, the corresponding tube waves overlap, which poses a particular challenge for the inversion process. Before the data are supplied to the inversion algorithm, we separated the P-wave from the tube waves, applied a move-out correction to the P-wave and then calculated a mean trace. A time-gated version of this mean trace with a length of 10 ms then serves as the estimate of the source wavelet. For this problem with three fractures, we have 15 unknowns, which are specified in Table 1. Three unknowns are related 185 to the background rock. These are the bulk and shear moduli of the formation and a separate shear modulus used for the tube waves. As outlined above, we use separate shear moduli for the P-and for the tube waves as a first-order approximation to account for anisotropy, which was estimated to be approximately 10% at the considered site (Wenning et al., 2018). The next Table 1. Unknowns of the inverse problem and their prior ranges subdivided by horizontal lines into three groups. The first group from the top comprises the background medium parameters, the second group the fracture parameters, and the third group algorithmic "tuning" that the fracture inclination cannot be inferred with high confidence, we assume that the inclination is known from televiewer data. The remaining three unknowns are algorithmic "tuning" parameters without any physical meaning. The first parameter of this group is the source depth. While the actual source location is known, we estimate the source depth for a fictitious homogeneous background medium to accommodate possible variations of the background medium parameters above the section under consideration. If the background rock is indeed homogeneous, the estimated source depth will correspond to the true 195 source depth. The other two "tuning" parameters are used to emulate attenuation of the tube waves. As tube waves propagate along the borehole, they do not suffer from geometrical spreading as, for example, the P-wave does ( Figure 3). However, tube waves are attenuated due to inelastic effects or scattering. To account for this, we dampen the tube waves using an exponential function defined by a shift factor, which specifies when the damping starts, and an exponent, which specifies the damping rate.
We ran our algorithm three times to ensure that it successfully locates the posterior PDF and does not get stuck in a local 200 minimum. Each time, three parallel Markov chains were used to explore the parameter space. More chains would have allowed for a more comprehensive exploration of the solution space, but would also have required more computational resources. Three chains are in our experience sufficient to exhaustively explore a 15-dimensional solution space well, such that the posterior PDF is found in most of the runs. The development of the root mean square error (RMSE) is plotted in Figure 4 for each that not all the data can be explained by the proposed model. With the objective to force the algorithm to more extensively explore the posterior distribution, we fix the standard deviation of the data error at a relatively high value, which is larger than corresponding estimates obtained in previous inversion runs. Figure 4 shows that all runs converge to a stable RMSE-value, which, as the data error is fixed at a high value, is smaller than one. Before reaching a stable RMSE, the algorithm explores the 210 complete solution space in search of the posterior PDF. This is referred to as the burn-in phase. Subsequently, the algorithm is expected to have located the posterior PDF and to explore it in the course of the remaining iterations. Steps per chain burn-in phase exploration phase Figure 4. RMSE weighted by the standard deviation of the data error for the three inversion runs of the observed VSP data shown in Figure   3. As the estimate of the standard deviation of the data error is fixed at a high value, the RMSE drops below one. The vertical black line indicates the separation of the burn-in and the exploration phases, associated with the MCMC search of the parameter space.
In order to assess whether the Markov chains have converged sufficiently to allow for a reliable estimation of the posterior PDF, we calculate the so-called potential scale reduction factorR (Gelman and Rubin, 1992). Considering only the part of the Markov chains after burn-in,R compares the variance of the individual Markov chains with the overall variance of all the 215 chains together. Usually, convergence is considered to be reached ifR is smaller than 1.2 for all parameters. In this example, considering a burn-in phase of 30% of the complete chains, we getR < 2 for most parameters, but only approximately a third of the parameters reachR < 1.2. Consequently, the posterior PDF has not been fully explored. Therefore, we do not plot posterior PDFs for the inferred parameters. Instead, we show the development of the Markov chains as a function of iteration number.
Although proper convergence has not been achieved, the inferred models explain the data well. However, other models, not 220 sampled by the Markov chains, might explain the data equally well. Hence, longer chains would be necessary to ensure a comprehensive exploration of the posterior PDF.
The acceptance rate specifies how many of the tested models are accepted. A too high acceptance rate generally implies that only models in the immediate neighborhood of the current model are explored while a too low acceptance rate means that computational resources are wasted by testing unrealistic models. Ideally, the acceptance rate ranges between 10 and 30%. In 225 our case, it lies between 10 and 20% for runs one and two and around 5% for run three. The most interesting inferred parameters are the apertures and compliances of the fractures, and to a lesser extent, the background rock properties. The development of these unknowns as a function of the number of iterations is plotted for all three runs in Figure 5. For the aperture of the first fracture (Figure 5a), the algorithm either finds a very large value of 10 mm (run 1) or a rather small one of less than 0.5 mm (runs 2 and 3). Interestingly, the opposite is the case for the second fracture 230 (Figure 5b). Here, run 1 suggests a small fracture aperture and runs 2 and 3 a large one. As mentioned earlier, the first two fractures are very close together, at 23.5 and 23.9 m depth, respectively. Hence, the corresponding tube waves overlap. The algorithm, thus, finds that one fracture must have a much larger aperture than the other, but it cannot determine which one is which. This leads to a bimodal posterior PDF featuring two equally probable modes. The estimated compliance values for these two fractures behave similarly (Figures 5d and 5e), although the difference between the runs is smaller.

235
The vertical axis of the plots in Figure 5 represents the prior range. In the cases where the first or the second fracture is found to have a large aperture, the inferred value is actually located at the upper limit of the prior range. This means, that the algorithm would propose even larger values if it were allowed to do so. We have not extended the prior range, because (1) even larger fracture apertures seem unrealistic and (2) the models found within this prior range are able to explain the data well.
The posterior PDF for the estimates of the aperture of the third fracture is unimodal (Figure 5c). At the location of the third 240 fracture, televiewer data (Krietsch et al., 2018) also indicate the presence of a larger shear zone. As we were not sure if the observed tube wave stems from the shear zone or the fracture, we extended the prior range of the aperture for this fracture by one order-of-magnitude to be able to accommodate the complete shear zone. However, all three runs suggest a small aperture of less than 1 mm, which clearly indicates that the tube wave is generated by the fracture and not by the shear zone.
For the bulk and shear modulus of the background (Figures 5g and 5h), we observe a similar behavior as for the fracture 245 apertures of the first and the second fracture: If the bulk modulus is large, then the shear modulus is small and vice versa. Both parameters are constrained by two observables: (1) the P-wave velocity by the moveout of the P-wave and (2) the transmission coefficient by the amplitude loss of the P-wave across fractures. However, these two observables are insufficient to constrain the background moduli adequately, thus leaving some degree of ambivalence in the final estimates. Conversely, the shear modulus used for the calculation of the tube-wave velocity is well constrained (Figure 5i), because there is no trade-off with other 250 parameters.
As the RMSE in Figure 4 is the same for all runs, the two modes of the posterior PDF identified by the algorithm explain the data equally well. To further illustrate this, we compare in Figure 6 synthetic data based on the inversion results presented in Figure 5 with the observed data. We generate the synthetic data using the last model of the third Markov chain of run 1 (blue in Figure 6a), in which the first fracture is inferred as having a large aperture, and of run 2 (red in Figure 6b), in which the second fracture has a large aperture. The observed data are plotted in black. Although we use a semi-analytic forward solver, which is inherently subject to a number of rather stringent assumptions, such as, for example, a homogeneous background medium, both synthetic datasets fit the observed data remarkably well.

Discussion
Observed data Simulated data: Run 1, chain 3 Observed data Simulated data: Run 2, chain 3 a) b) Figure 6. Comparison between simulated (colored) and observed (black) data: (a) run 1 and (b) run 2 study apertures of 6.4, 1.7 and 0.0 mm. These are the fracture apertures at the borehole wall, which are not identical to the hydraulic fracture apertures inferred in this study. While the former represents the actual aperture at the interface between the borehole and the fracture, is the latter an average of the hydraulic aperture over the rock volume in the vicinity of the borehole sampled by the VSP data. In spite of these differences, the televiewer data can, for example, help identify the correct mode in case of multimodal posterior distributions. In our case, the televiewer data indicate that the first fracture has a larger aperture 265 than the second one, confirming that the modal aperture distribution identified by run 1 is realistic. However, run 1 infers for the second fracture a much smaller aperture than indicated by the televiewer data. This indicates that, although the fracture has according to the televiewer an aperture in excess of 1 mm at the borehole wall, it is likely to be much thinner away from the borehole. The aperture of the third fracture is smaller than the vertical resolution of the optical televiewer of 0.21 mm. Similarly, we also obtain a very small fracture aperture, that is estimated by our algorithm to be below 1 mm. Concerning the fracture 270 compliances, we can compare our results with those of Barbosa et al. (2019), who present corresponding estimates for the same borehole section based on full-waveform sonic log data. They estimated fracture compliances which are approximately one order-of-magnitude higher than our results (9.9 · 10 −13 m/Pa). Potential reasons for this difference might be that the fullwaveform sonic data were measured at significantly higher frequencies (∼ 20 kHz) than our VSP data and that the fracture compliances tend to be frequency-dependent (e.g., Pyrak-Nolte, 1992;Nakagawa, 2013). Another difference between the two 275 studies is the incidence angle. While Barbosa et al. (2019) assume normal incidence of the P-wave on the fractures, this study accounts for the dip angle of the fractures derived from televiewer data, which ranges from 62 • to 77 • with regard to the horizontal.
A bit puzzling is the remarkably low estimate of the tube-wave shear-modulus of only about 6 GPa ( Figure 5i). This parameter is very well constrained, as it is the only free parameter in equation 5, which may, however, be too simplistic for the 280 following three reasons: (1) Equation 5 is derived in the low-frequency regime and its validity for higher frequencies is limited.
(2) Attenuation of tube waves, as for example through scattering on the borehole tool or inside the damaged zone surrounding the borehole, was not accounted for when estimating the tube-wave shear modulus.
(3) Anisotropy is not taken into account completely. Thus, while the resulting tube-wave velocity is correct, as can be seen by the excellent fit between the observed and synthetic data, the corresponding shear modulus appears to be underestimated in order to correct for physical effects neglected 285 in equation 5. Incorporating attenuation into the tube-wave velocity equation can be done by implementing equation 5-17 of White (1983) including the impedance of the borehole wall and accomodating anisotropy can be done by one of the methods presented by Karpfinger et al. (2012). This, however, is beyond the scope of the present study.
From an inversion perspective, the most interesting point of these results is that two modes of the posterior PDF were identified. This showed that having the first fracture with a large aperture, while the second fracture is thin, is similarly probable 290 as the opposite scenario. Note that a deterministic approach would have provided only one result without any indication that there is another mode that can explain the data equally well, whereas our Bayesian approach clearly supplied us with both options. This nicely demonstrates the value of stochastic inversion approaches.
A downside of our Bayesian approach is its enormous computational cost. Most of it is spent in the forward steps to simulate VSP data for the proposed model. We have optimized the forward simulation by parallelizing over frequencies. Still, one 295 inversion run with three parallel Markov chains and 60'000 forward steps per chain took approximately 14 days to complete using one node (48 AMD Opteron 6174 processors at 2.2 GHz) of our cluster. However, the inversion would run three times faster if each of the three Markov chains were run on a different node. We did not do this due to limited availability of resources.
In any case, we argue that the computation time is well spent, since the results obtained are much more comprehensive than those that would be obtained through a deterministic inversion, as they allow, as explained above, to discover multiple modes 300 of the posterior PDF. Furthermore, stochastic inversion approaches do not really depend on the starting model. This is in stark contrast to deterministic full-waveform inversion approaches, which require starting models whose forward response deviates from the forward response of the true model less than half a wavelength (Virieux and Operto, 2009).
For the real-data example, we have decided not to estimate the source wavelet during the inversion process, although the corresponding algorithm was developed and successfully applied for synthetic test cases as demonstrated in the first results-305 section. The reason is, that the source wavelet of the observed data includes extensive reverberations and is, thus, extremely long and complicated. Estimating it as part of the inversion procedure would have required to more than double the amount of unknowns, which would have rendered the problem unnecessarily complex.
An important limitation of our forward model, and indeed of virtually all fracture-based tube wave models, is that fracture aperture and compliance are correlated. This means that the inversion algorithm tends to compensate an overestimation of the 310 fracture aperture by underestimating the fracture compliance. Therefore, we observe that a large fracture aperture for the first fracture is accompanied by a relatively small fracture compliance (Figures 5a and 5d). This is supposed to be mitigated in our approach, because the estimate of the fracture compliance is not only constrained by the tube-wave amplitude, but also by the reduction of the P-wave amplitude when a fracture is crossed (Schoenberg, 1980). However, the transmission coefficients calculated for the estimated parameters are very close to 1, and hence, the effect of this constraint is relatively weak. As the 315 Markov chains are not oscillating all over the prior range, and as the obtained values are reasonable, we can conclude that this compensation is rather limited.
Inspecting the difference between the observed and the forward modeled data shows that the largest discrepancies are found at the fracture locations. This indicates, that the transmission loss of the P-wave across fractures may not be reproduced properly in the synthetic data. However, as this affects only the P-wave around the fracture locations, the impact on the RMSE 320 are limited. A possible way to improve this issue might be to define a weighting function that peaks at the fracture locations to force the algorithm to obtain a better data fit at these locations, and thus, find a more accurate transmission coefficient. The downside of this, however, is that the weights are new "tuning" parameters that need to be adjusted through a time-consuming process, which was not feasible to accomplish in the scope of this study.
Limitations of our implementation of the forward operator are its inability to account for scatterers, impedance contrasts 325 related to lithological changes, and borehole enlargements. If corresponding effects are present in the data, they might need to be filtered out prior to inversion. Similarly, changes in the P-wave velocity are not taken into account. If these are present, the data needs to be cut into smaller pieces with constant P-wave velocity. Changes in P-wave velocity above the considered borehole section are taken into account by virtually shifting the source depth. The algorithm is also not able to take S-waves and corresponding tube waves into account. In our dataset, events of this kind were indeed present and needed to be muted 330 before applying our inversion algorithm to it.

Conclusions
We have developed a Bayesian MCMC full-waveform inversion algorithm based on a semi-analytical forward solver to simultaneously infer the aperture and compliance of individual fractures from corresponding tube-wave data. We mitigate the correlation between fracture aperture and compliance by constraining the fracture compliance by two independent observables:

335
(1) the tube-wave amplitude relative to the P-wave amplitude and (2) the amplitude loss of the P-wave across a fracture. The algorithm was applied to a field dataset acquired in crystalline rock at the underground Grimsel Test Site in Switzerland. The subsection of the VSP dataset considered contained three fractures, of which two are very close together. The algorithm identified two equally probable modes in the posterior PDF: Either the first fracture features a large aperture and the second fracture a small one or vice versa. In other words, from the information provided, the algorithm can determine that one fracture is 340 larger than the other, but it cannot determine which one is thick and which one is thin. The identification of these two modes clearly illustrates a major advantage of stochastic inversion algorithms as compared to their deterministic counterparts. The latter would not have identified these two modes and would have provided just one of the two possible solutions. Our case study also shows that in a complex geological environment with multiple, closely spaced fractures, the hydraulic apertures of individual fractures cannot be determined. However, the method can still provide an effective fracture aperture distribution of 345 a package of fractures. The inferred apertures in our example are consistent with televiewer data and the inferred compliances are roughly in the same range as those derived from sonic logs at the same site. The data fit is remarkably good, especially when considering the semi-analytical nature of the forward solver and the inherent assumptions it relies on, as well as, the rather complex character of the observed hydrophone VSP data.