T H E SPATIAL FILTERING F R O M C O R E - S C A L E C O N D U C T I V I T Y TO T H E C O N D U C T I V I T Y M E A S U R E D B Y A S L U G TEST By Bin Wang B. Eng., Hydrogeology and Engineering Geology Xiangtan Mining Institute M. Sc., Hydrogeology and Engineering Geology Chinese Academy of Geological Sciences A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF M A S T E R OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES GEOLOGICAL SCIENCES We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA August 1995 © Bin Wang, 1995 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Geological Sciences The University of British Columbia 6339 Stores Road Vancouver, Canada V6T 1Z4 Date: Abstract In the thesis, I investigate the measurement process of the slug test, an important and efficient engineering technique employed to probe aquifer properties. The objective is to determine the scale at which a slug test 'measures' hydraulic conductivity. To determine the scale, I conceptualize slug test measurement process in a system and filtering frame-work. I quantify, through numeric experiments, the spatial filtering effect of the slug test upon the core-scale hydraulic conductivity field. I examine the role of specific storage on this filtering effect. I develop two approaches to identify the system filter that epitomizes the filtering system — measurement process: spectral analysis with Wiener filtering and numerical perturbation. In heterogeneous media, system linearity is analyzed and the inverse es-timation of the system output is evaluated. With the spectral approach, I optimally estimate the system filter in the presence of measurement noise. In homogeneous media where the spectral approach cannot be used, I apply the numerical perturbation ap-proach. The spectral analysis and numerical perturbation are two conceptually different but complementary approaches. I analyze the nonparametric form of the system filter in heterogeneous media of var-ious spatial variability and specific storage values. I determine the empirical parametric expression of equivalent filter width versus specific storage in homogeneous media. It is found that: 1. In both heterogeneous and homogeneous media, the equivalent filter width increases linearly with respect to log 2. In homogeneous media, the filter amplitude decreases away from the wellbore by an approximate ^ law if Ss < 10 - 4[L - 1]; ii 3. A preliminary analysis of the role of heterogeneity characteristics upon the slug-test filtering was undertaken. The slug-test filtering was not strongly influenced by heterogeneity characteristics for log-conductivity variances less than 0.7. At variances less than 0.3, anisotropy in the conductivity field caused anisotropy in the slug-test filtering. These effects were not quantified in this thesis. In the end, I discuss the possible engineering and theoretical applications of the spatial filtering approach and make recommendations for future work. iii Table of Contents Abstract ii List of Tables ix List of Figures x Preface xvii i Acknowledgment xix Dedication xx 1 I N T R O D U C T I O N 1 1.1 Scale and averaging 1 1.1.1 An opening example 1 1.1.2 Scale and averaging: A filtering perspective 2 1.1.3 Parameter upscaling: A filtering interpretation 11 1.1.3.1 Parameter measurement: a special case of upscaling . . . 12 1.2 The measurement filter: slug test 15 1.2.1 Review of slug test literature 15 1.2.2 The slug test filter 17 1.3 Thesis overview 18 2 S I M U L A T I O N O F S L U G T E S T 19 2.1 The boundary value problem of a slug test 21 2.2 2D FDM in rectangular coordinates 23 iv 2.2.1 Numerical approximation and discretization errors 24 2.2.1.A Wellbore approximation 24 2.2.l.B Spatial discretization 26 2.2.1.C Time discretization 31 2.2.1.D The hypothetical remote boundary 31 2.3 2D FDM in radial coordinates 33 2.4 ID FDM in radial coordinates 34 3 S Y S T E M A N D F I L T E R I N G 38 3.1 The concepts of system and filtering . 38 3.2 The slug test conductivity filtering system 39 3.2.1 Physical definition 39 3.2.2 Mathematical representation 44 3.3 System inputs: core-scale fields 45 3.4 Properties of System filtering 46 3.4.1 Stability and nonuniqueness 46 3.4.2 Space invariant property 46 3.4.3 Linearity 48 3.5 System outputs: inverse estimation of K and Ss 50 3.5.1 Objective function versus estimation variables 53 3.5.2 Comparison of the estimates of different inversion methods . . . . 57 3.5.3 Interpretation of the estimates of different inversion methods . . . 62 3.5.3.A Extreme values of conductivity estimated with true S3 . 63 3.5.3.B Effects of spatial conductivity variability with different specific storages on the inverse estimates 68 v 3.5.3.C Influence of low and high conductivity at different loca-tions from the wellbore on the inverse estimates 77 4 S Y S T E M I D E N T I F I C A T I O N : S P E C T R A L A N A L Y S I S 83 4.1 Nonparametric filter function identification: Wiener filtering 83 4.2 Spectrum estimation 85 4.2.1 Bias and dilemma of spectrum estimation 85 4.2.2 Periodogram and data tapering 89 4.2.3 Multitaper spectral analysis 90 4.2.3.A Evaluation: variance and bias 92 4.2.3.B Verification: hypothetical case 95 4.3 Testing the Wiener filtering method 99 5 S Y S T E M I D E N T I F I C A T I O N : P E R T U R B A T I O N M E T H O D 102 5.1 Introduction 102 5.2 Numerical perturbation formula 104 5.2.1 The inverse estimation procedures of perturbed conductivity . . . 109 5.3 The area and the amount of perturbation 115 5.4 Filter functions of different specific storages 119 5.4.1 Fixed perturbation magnitude 120 5.4.2 Ss-dependent perturbation 121 5.4.3 Parametric representations 131 6 F I L T E R F U N C T I O N 134 6.1 Filter function for different specific storages 134 6.1.1 Slug tests in isotropic core-scale fields: different field sizes, different specific storages 135 vi 6.1.2 Comparisons of Kcore and Kalua spectra 142 6.1.3 Comparisons of filter spectra, filter functions 153 6.2 Filter functions in conductivity fields of different spatial variabilities . . . 168 6.2.1 Variance, variogram model and stochastic realization 168 6.2.1. A Multiple realizations of isotropic fields 168 6.2.2 Core-scale conductivity of different correlation length 172 6.2.2. A Isotropic, different correlation length 172 6.2.2.B Anisotropic, different correlation length ratios 175 6.2.2.C Lense field 183 7 C O N C L U S I O N S 188 7.1 The spatial filtering from core-scale conductivity to the conductivity mea-sured by a slug test 188 7.2 Possible engineering application 190 7.2.1 The Q-problem 190 7.2.2 The 2,-problem 191 7.3 Theoretical significance 192 7.4 Future work 192 Bibliography 194 Appendices 201 A E F F E C T I V E A N D E Q U I V A L E N T P A R A M E T E R S 201 A.I Effective parameter: a limiting case of upscaling 202 A.II Equivalent parameter: a general case of upscaling 207 vii B F I N I T E D I F F E R E N C E E Q U A T I O N S O F S L U G T E S T 209 B.I 2D FDM in rectangular coordinates 209 B.II 2D FDM in radial coordinates 214 B.III1D FDM in radial coordinates 217 viii List of Tables 3.1 Linearity analysis of K, \gK and 49 3.2 The maximum and minimum estimates, skewness K*, AK and ASS. . . 63 3.3 K", AK ~ ASa of tests a, b, c, d. 76 5.4 The G(l + (), p of the filter expression G(rD) = G(l + 0 T - - ^ 131 6.5 Slug tests with different specific storages 135 6.6 Slug tests in isotropic core-scale fields from multiple realizations 169 6.7 Slug tests in anisotropic conductivity fields 175 6.8 Statistic parameters for two lense fields 183 6.9 Slug tests in two lense fields 183 ix List of Figures 1.1 Scale and averaging: real world, engineering world, model, and decision analysis 4 1.2 The A(scale) ~ a(correlation range) of continuously filtered real world(A), measurements(C), model(E), decision analysis(F), and filtering processes(B, D) 9 2.3 Domain discretization of slug tests. Graph A is a rectangular discretiza-tion in heterogeneous conductivity field; B is a radial discretization in homogeneous field with a small region of inhomogeneity(only half shown); C is a radially symmetric discretization in homogeneous field 20 2.4 A slug test in a confined aquifer 21 2.5 A: Discretizing wellbore into a number of fine blocks. B: A round well is approximated by a square block 25 2.6 Head distribution at initial time and head response at wellbore 26 2.7 Wellbore response comparison: numerical solution against analytical solu-tion, a: space step equal a wellbore diameter; b: space step equal ^ of a wellbore diameter. All other parameters of a and b are identical 27 2.8 Wellbore response comparison: numerical solution against analytical so-lution. The solid curves are analytical solutions; The dashed curves are simulated results. Specific storage ranges from 10~4(upper right) to 10 - 1 (lower left). Their space steps are | , | , of a wellbore diameter respec-tively 28 2.9 Effective screen radii in terms of wellbore diameter versus logarithmic specific storages 29 x 2.10 Effective screen radii in terms of wellbore diameter(Ao;u;) are a function of time. The smaller S3, the flatter the re(t), the faster re(t) reach 0.564Axu'. (The last time step is, t55 = ^ = ""P-™**")') 3 0 2.11 Simulated wellbore response(dashed lines) versus analytical solution(solid lines). The graph only shows the curves of S„ ranging from IO - 5 (upper right) to 10 - 1 (lower left). For S3 < 10 - 5, the fitting improves 35 2.12 Simulated wellbore response (dashed lines) versus analytical solution(solid lines). From lower left to upper right, S3 decreases from 10 _ 1 to 10 - 1 0 . For all simulated curves, 8 = = 10 37 3.13 Slug test filtering system receives core-scale conductivity field and pro-duces measured-conductivity field on a larger scale 40 3.14 Slug test scenario 41 3.15 Spatial conductivity filtering is a transient process. The top three show the core-scale field under filtering at different time t\, t2 and t^, where tz > t2 > t\. The bottom are the corresponding filter functions 43 3.16 Scaling property Y®: (Yc -> Ym) ==• (a Yc -» a Ym), where a\c = 0.01. . 51 3.17 Additivity y®: (Ycl Yml, Yc2 ^ Ym2) (Yel + Yc2 -• Yml + Ym2). . . 52 3.18 Objective function versus logarithm conductivity and specific storage: e( lg K, lg S.) versus (lg K, lg S.). Top graph, e(lg Ktrue, lg Sf u e ) = e(2, -1) = 0; bottom graph, e(\gKtrue, lgS' r u e) = e(2, -10) = 0. Contour level: 0.4, for e(lgK, \gSs) > 0.4; 0.0025, for e(\gK, lgSs) < 0.2 55 3.19 Comparison of histograms of the IgK estimated with true lg^s and the ones of the \gK concurrently estimated with lg Ss 59 3.20 Cross-plots of lg K concurrently estimated with lg Ss versus estimated \gS3. The dashed lines represent the true lg5 a 60 xi 3.21 Comparison of lg K estimated with true lg Sa with lg K concurrently esti-mated with lg Ss- The inclined lines are of unit slope 61 3.22 Cross-plots of [(lgi»f estimated with true \gS„) - (IgK estimated concur-rently with lg •!?«)] versus [(estimated \gS3) - (true lgS*)] 62 3.23 The solid lines are the simulated response curves. The dashed lines are the fitting response curves. Left column represents the fit by varying K only. Right column is the fit by varying both K and Ss 67 3.24 Conductivity variation functions and skewness K** of tests a, b, c, d. . . 71 3.25 AK ~ ASa of conductivity variation functions a, b, c, d. 72 3.26 Fitting and response curves of K>, KK and K° of test d 75 3.27 Perturbation set-up, and the repsonse curves of low/high K perturtbations near wellbore and distant from wellbore 80 3.28 K± ~ r D , ~ rD, where rD = r/rw, of inversion scheme a and b. True Ss = 0.1 81 3.29 AK± ~ rD, ASS ~ rD, where rD = r/rw, of inversion scheme a and b. True Ss = 0.1 82 4.30 Filter identification by an iterative Wiener filtering procedure 86 4.31 Spectrum estimation procedure by spectral method 87 4.32 Variance, leakage and total spectral energy with respect to spectral window width, P = NW, and the number of prolate tapers, K, of the multitaper (left column) and adaptive multitaper methods (right column).A2B2 is the leakage for both methods 94 4.33 Comparisons of multitaper spectra, periodogram and true spectra. The primitive field is filtered by a Gaussian function G(x, y)=v £ A exp[—6(|j + ^)], where filter width Xx = Xy= 8, 16, 32Ax 98 xii 4.34 Left column: comparisons of the estimated filter function with the hy-2 2 ^ ( 2 I 2 ) pothetical filter function G(x,y) = v £ x^ e *x A» ; right col-umn: comparisons of the estimated transfer function with the hypothet-ical transfer function G(fx,fy) = e ~~ « ^ ^ A* + A» ); where filter width A x = A y = 8,16,32Ax respectively. . . 101 5.35 The radial mesh for the simulation of perturbation experiments 106 5.36 The IgK inverted by different methods and the associated filter functions. 113 5.37 Spatial filter is a function of time. A shows the filter function at locations Pi, P2, P3 at time [0,3], 8 = 50. B shows the filter function at time h, ^2, <3 in dimensionless space rD, rD € (1, Re] 114 5.38 Perturbation with different dk. The \gK is estimated with lgS", equal to true value 117 5.39 Perturbation with different dk. The IgK is concurrently estimated with \gSa. 118 5.40 The spatial filter for different specific storages, Sa = 10 - n , n = 1,2,..., 10. For all perturbations, dk = 0.5. A shows the filter function versus dimen-sionless distance in log space. B is the volume of filter function integrated Re over dimensionless distance rD, VG(V d) = / 2 w G(rD) rD drD 123 5.41 A is the equivalent filter width defined as Ae = cffi^j • B is the influence radius Re defined as -> 0 124 5.42 Sa = 10~x. A shows that G{rD) ~ rD converges to the dk = 0.01 curve. If dk is too small(t4 < 0.01), solutions become instable. B is the correspond-ing VG(rD) ~ lg(^j) curves 125 xiii 5.43 Ss = IO - 3 . A shows that G(rD) ~ rD converges to the dk = 0.05 curve. If dk is too small(<4 < 0.05), solutions become instable. B is the correspond-ing VG(rD) ~ lg(^j) curves 126 5.44 Ss = 10 - 5. A shows that G(rD) ~ rD converges to dk = 0.1 curve. B is the corresponding VG{TD) ~ lg(^j) curves 127 5.45 Sa = 10 - 7. A shows that G(rD) ~ rD converges to dk = 0.3 curve. B is the corresponding vci(rD) ~ lg(^j) curves 128 5.46 A shows the Va(rD) ~ lgr D curves calculated with the optimal dk for Ss = 10 _ n, n = 1,2,..., 10. B is a comparison of the Ae ~ lg(^-) of A with those estimated with dk = 0.5 for all Ss values. C is the corresponding k[G(rD)] ~ \grD 130 5.47 A: Fitting G(l + 0 ~ ^ by G(l + <t) = 0.25s025 + 0.25J-76 + 0.01753; B: Fitting p ~ ^ by p = 3.0y/Ss + 2.0 132 6.48 The grey-scales and histograms of core-scale and measured fields with different specific storages. The core-scale field has an isotropic exponential variogram. Field size 64 x 64 blocks 137 6.49 The variograms of core-scale and measured fields of size 64 x 64 138 6.50 The variograms of core-scale and measured fields. Field size 128 x 128. . 139 6.51 The variograms of core-scale and measured fields. Field size 256 x 256. . 140 6.52 The mean and variance of core-scale and measured fields with different specific storages. Field size 64 x 64, 128 x 128, 256 x 256 141 6.53 Spectra of the core-scale fields and the measured fields with different spe-cific storages. The core-scale field has an isotropic exponential variogram with a\K = 0.25, A* = Ay = 4Ax w 146 xiv 6.54 Comparisons of the spectra of core-scale fields and measured fields of size 64 x 64, 128 x 128, 256 x 256 blocks 147 6.55 Spectra of the core-scale fields and the hypothetically filtered fields with different filter widths. The core-scale fields are generated with a sequencial Gaussian method of GSLIB, with an exponential variogram, tr^^- = 0.25, \ x = \ y = 4Ax w 148 6.56 Comparisons of the spectra of the core-scale field generated with GSLIB and the hypothetically filtered fields of size 64 x 64, 128 x 128, 256 x 256 blocks 149 6.57 Spectra of the core-scale fields and the hypothetically filtered fields with different filter widths. The core-scale field is generated with a spectral method with an exponential variogram, o~fsK = 0.25, Xx = \y = AAxw. . . 150 6.58 Comparisons of the spectra of the spectrally generated core-scale and the hypothetically filtered fields of size 64 x 64, 128 x 128, 256 x 256 blocks. 151 6.59 The Rayleigh freuencies and the target spectra of core-scale and filtered fields 152 6.60 Filter spectra of the slug tests in the conductivity fields of different specific storages 158 6.61 Comparisons of the filter spectra in fields of size 64 x 64, 128 x 128, 256 x 256 blocks 159 6.62 Averaged filter functions for different specific storges. . 161 6.63 Filter functions for different specific storges from spectral method. En-larged by 104. Contour level = 200 162 6.64 Filter functions for different specific storges from perturbation method. Enlarged by 104. Contour level of inner zone = 200; contour level of outer zone = 0.2. . 163 xv 6.65 Scatter plots of measured fields versus estimated fields from Wiener filter-ing. Field size 256 x 256. The inclined lines are of unit slope 164 6.66 Scatter plots of measured fields versus estimated fields from Wiener filter-ing. Field size 128 x 128. The inclined lines are of unit slope 165 6.67 Scatter plots of measured fields versus estimated fields from Wiener filter-ing. Field size 64 x 64. The inclined lines are of unit slope 166 6.68 Comparison of the equivalent filter width versus specific storage (Ae ~ lg(-^j)) from the spectral analysis and the perturbation method 167 6.69 Filter functions estimated in multiple core-scale conductivity realizations, with different variances and variogram models. The last graph is a com-parison of the averaged filter functions shown in the first three graphs. . 171 6.70 Comparions of filter functions in isotropic fields of different correlation lengths, A x = Aj, = 3, 6, 12Axw 173 6.71 Scatter plots of the measured fields versus the estimated fields from Wiener filtering. The core-scale fields have isotropic vaiograms of correlation lengths A = 3, 6, 12Axw 174 6.72 The grey-scales and histograms of the core-scale and the measured field of different correlation length ratios. The fields shown are 64 x 64 blocks each 178 6.73 Graph A: reduced variance(«7yc — o-Ym) versus correlation lenght ratio Graph B: the mean of measured fields fiYm versus the mean of core-scale fields HYC 179 6.74 The filter functions on the center lines in the X and Y directions in anisotropic conductivity fields 180 6.75 Filter functions of different anisotropy ratios. Enlarged by 104 181 xvi 6.76 Scatter plots of the measured fields versus the estimated fields from Wiener filtering. Anisotropy ratio of the core-scale fields: \ x / \ y = 1, 3, 5, 30Axw. The inclined lines represent unit slope 182 6.77 The grey-scales and histograms of the core-scale and the measured field with lense heterogeneties. Size of Lense Field one = 128 x 128; size of Lense Field two = 64 x 64 185 6.78 Poorly resolved filter functions in two lense fields. Enlarged by 104. . . . 186 6.79 Scatter plots of the measured fields versus the estimated fields from Wiener filtering. The inclined lines are of unit slope 187 B.I.80 Block centered finite difference in Cartesian coordinates 209 B.I.81 Wellbore specific stoarge equals to 1 212 B.II.82 Mesh centered finite difference in radial coordinates 215 B.II.83 Finite difference of the inner boundary. 216 xvii Preface The thesis subject is the spatial filtering from core-scale conductivity to the conductiv-ity measured by a slug test. In Chapter 1, we introduce a spatial filtering approach to interpret a slug test measurement. Some relevant concepts are discussed from a filtering perspective. In Chapter 2, we present numerical methods to simulate slug tests in het-erogeneous and homogeneous media. The conceptualization of the slug test in a system and filtering framework is given in Chapter 3, where system filtering properties and the inverse estimation of system outputs are discussed in detail. In Chapter 4, a spectral analysis through Wiener filtering is developed to identify spatial filter in heterogeneous media. As a result of computation load constraints, the specific storage value examined in the heterogeneous conductivity filtering is limited to Ss = 10 _ n[L _ 1], n = 1, . . . , 5. The effects of other Ss with values of n = 6, . . . , 10 on the conductivity filtering are investigated in homogeneous media where the spectral analysis method becomes invalid. A numerical perturbation method is developed, which is the topic of Chapter 5. In Chap-ter 5, empirically determined expressions of influence radius and equivalent filter width versus specific storage are presented. Consequently, Chapters 4 and 5 contain conceptu-ally different but complementary approaches. In Chapter 6, the conductivity filtering in heterogeneous media with different specific storages and various spatial variabilities are explored. In Chapter 7, the major results are summarized in conclusions and future work is recommended, also the possible engineering application and theoretical significance of the spatial filtering approach are discussed. xviii Acknowledgment First of all I wish to thank my advisor, Roger Beckie. Without his help this thesis would not have been written. His help always extends to more than research, and I appreciate his various discussions on culture and history around the world. I would like to thank my committee member Les Smith. His inspiration on the practical applications of my thesis results is of great value. My committee member Tad Ulrych deserves my heart-thanks. His thorough review and comments on the system and filtering and on the spectral method leads my thesis to a more robust version. I am indebted to Greg Dipple for his kindness being my committee chair. I would like to acknowledge my fellow students in the Hydrogeology group. I also wish to thank my former advisor, Xibin Zhu, and my colleagues at the Chinese Academy of Geological Sciences. It was their help that provided me knowledge on some advanced topics in hydrogeology. Support from my family is of great help for the completion of the thesis. I am especially grateful to my parents whom I can always count on. My gratitude to my grandmother who passed away during my study abroad will be cherished deep in heart and last forever. My sisters, Lichun, Liyuan, Lixin and Lifeng receive special thanks for helping Mom and Dad all the time. If there were one thought that I could best express by words in this thesis, I would devote it to colleague and my wife, Yanhong; that thought would embrace my whole-hearted appreciation for all that she has suffered through, and my sincere thankfulness for all that she has done. xix Dedication This thesis is dedicated to the memory of my grandmother XIYING HOU (17 May 1915 — 4 December 1993) xx Chapter 1 I N T R O D U C T I O N In this chapter, we introduce the slug test measurement from a spatial filtering perspec-tive. We first discuss some relevant concepts: scale, averaging, and parameter upscaling, under a filtering framework. The literature review of the work on parameter upscaling is relegated to Appendix A. The more relevant work on the measurement — pumping test is reviewed in Section 1.1.3.1, with focus on the influence of the aquifer heterogeneity on the measured conductivity. In Section 1.2 we review the literature on the slug test and we introduce a spatial filtering approach to interpret the slug test measurement. A brief thesis overview is given in Section 1.3. 1.1 Scale and averaging In this section we discuss the scale and averaging in the real world and the engineering world from a filtering perspective, and we present a filtering interpretation to parameter upscaling. 1.1.1 A n opening example Fig. 1.1.1 shows a common scale problem. The shadowed zones schematically represent the sampling zones of well test, core test and tracer test, which all are the basic techniques utilized to probe aquifer properties. The grids represent an engineering model which is to be employed to interpret real-world processes and predict field problems. Now we see that different engineering sampling techniques measure aquifer properties on different 1 Chapter 1. INTRODUCTION 2 scales. Then the scale problem is: How can one util ize the measurement data in the engineering model? To solve the problem, one needs to understand the scale relationship between different measurements and the scale relationship between measurement and real world processes. In the thesis, we develop a spatial filtering approach that quantitatively characterizes the scale relationship from one scale to another. In the following, we first discuss the concepts of system and filtering, then we return to the scale issue. Engineering Model different aquifer tests sample conductivity on different scales 1.1.2 Scale and averaging: A filtering perspective A system can be viewed as any process that results in the transformation of signals. Thus a system has an input signal and an output signal which is related to the input signal through the system transformation [Oppenheim et al., 1983, p35]. If the trans-formation changes the relative amplitudes of the frequency components in a signal or perhaps eliminate some frequency components entirely, the transformation is referred to as filtering [Oppenheim et al., 1983, p397]. If the system is linear and space invariant, Chapter 1. INTRODUCTION 3 the filtering process can be represented in a convolution form,1 oo Y(x) = J H(x-x') X(x') dx', (1.1) —oo where X and Y represent the input signal and the output signal respectively, and H(x) is called filter function, its Fourier transform is called transfer function. The so-called system identification is to find the filter or transfer function, e.g., the impulse response of a linear time/space invariant system.2 The filter function epitomizes the system, hence characterizes the relationship between input and output. Mathematically, the convolution represents a summation or superposition. Physically, the convolution signifies an averaging process, e.g., the measured conductivity [Y(x)] by well test is somehow averaged by [#(«)] from the core-scale conductivities [X(x)]. In the thesis averaging is congruent in meaning to filtering. The real world exhibits variabilities on infinite scales, and it is frequently perceived to be nonergodic/nonlinear. The process in the real world can never be exactly reproduced in its counterpart, the engineering world, where scientific knowledge holds valid only on finite scales and only limited scale ranges are of practical concern. The engineering world is often assumed to be linear, stationary, ergodic and Gaussian [Newland, 1993, p80]. These assumptions are riot true in cases where sufficient measurements are not available and the measurements are obtained on small scales. However, if more measurements are performed and averaged to larger scales, we are heading for the engineering assumption. It is averaging process that bridges engineering knowledge with reality. Averaging is always linked with the issue of scale, and averaging itself is a transformation of scale — scaling, usually upscaling. The scale and averaging are shown in Fig. 1.1 and discussed in the following. ^ o t e that the formulas in the thesis are mixture of 2D and 3D, where a single argument appears, it should be considered a vector coordinate. The slug test results in the rest chapters are 2 dimensional. 2 " / " denotes "and/or". Chapter 1. INTRODUCTION 4 real world infinite scale natural u engineering \measurement estimation engineering world finite scale model interpretation resolving scale decision analysis decision-variable scale Figure 1.1: Scale and averaging: real world, engineering world, model, and decision analysis. There are two possible broad groups of averaging operating in the real world and engineering world: natural averaging and engineering averaging. Natural averaging is postulated to be the real-world averaging that involves scale transformation. For example, a regional-scale fluvial fan is a result of the transformation of small-scale sediments(e.g., igneous fragments, metamorphic clasts etc.), the averaging processes may involve water process, temperature and pressure and so on. Natural averaging may occur on geological time scales, thus it is hard to understand through experiments,3 but we can interpret it in a filtering framework. First, the aquatic sedimentary environment is represented by a system operating on different scales; second, the unsorted rock, i.e., the source, as input 3 Note that the natural averaging has meaning on spatial scale greater than that of representative elementary volume(REV) pertaining to the soil or porous media composed by mixtures of powders, where a fractal phenomenon seems have been discovered [Ghilardi et al. 1993, Katz and Thompson 1985]. These authors explain self-similar heterogeneity by a fractal behavior due to growth of crystals within the pore space and relate fractal dimension to morphogenetic kinetics. Here, the natural averaging operates on the sediments which are already-formed, apparently on larger scales than those of crystals. Chapter 1. INTRODUCTION 5 signals; last, the observed aquifer structures constitute system output. Consequently, the system filtering/averaging signifies the scale relationship between the source rocks and the aquifers. Natural averaging is beyond the scope of this thesis, we shall not discuss it further. The engineering averaging is of the most concern and it may be excited by various measurement techniques. A measurement process can be defined as, in a mathemati-cal notation, a mapping from the real-world process to an engineering framework, the notation has an identical basic form as (1.1), o o Y(x)p = J H(x - x')p X(x') dx' (1.2) —oo where Yp is a measure of X and Hp denotes the filter function associated with measure-ment process P [Baveye and Sposito, 1984; Cushman, 1984; Cushman, 1986]. We use a terminology, the equivalent filter width of the filter function, to represent the scale of measurement P, Ap, defined as [Bracewell, 1986] +°° / H(x)p dx X p = ~°° H(0)p ' ( L 3 ) that is, the area of the function divided by its central coordinate. Note that eq.(1.3) is a one dimensional expression.4 If time is included, eq.(1.2) extends to o o Y(x,t)p = J H(x- x',t - t')p X(x',t') dx'dt'. (1.4) —oo For example, a slug test 'measures' different aquifer volume at different times, the conductivity mapped from the core-scale(wellbore-scale) conductivities, by model of 4 Note that the equivalent filter width differs from filter width. If we consider filtering as a random signal process, filter width(A) has the same meaning as the correlation range(a) of the signal's variogram, and as the integal scale of its covariance function(r), i.e., a = r = X. In this thesis, we save the term 'filter width ' to describe filters, while 'correlation range' for variogram and 'integral scale' for covariance of random signals. Chapter 1. INTRODUCTION 6 Cooper et al. [1967], is a volume and time average. As noted by Cushman [1986], one can not 'measure' constitutive variables, e.g., conductivity, storativity, and one can only measure fundemental variables, e.g., hydraulic head, flux. Constitutive variables can be computed from the constitutive equations in a closure model, e.g., the model of Cooper et al. [1967].5 This averaging of (1.2) and (1.4) is namely spatial or spatial-time averaging, more specifically, we call it measurement averaging or filtering. In practice, all the variabilities can never be entirely observed, since this would require enormous measurements up to the sampling volume covering, for instance in mining engineering, the field of variability from lpm to lOOArm [Journel and Huijbregts, 1978, pl49]. Hence, measurements are only taken at some selected localities pertaining to some optimal sampling network with a risk-cost-benefit objective [Freeze et al., 1990; Bogardi et al., 1985]. The spatial structure of the measured variable inferred from measured data is then utilized in a variety of geostatistical, statistical or interpolation methods [Deutsch and Journel, 1992; Cressie, 1991; Journel and Huijbregts, 1978; Marsily, 1986], to estimate the variable state in unsampled region. This type of estimation bears an essence of averaging,6 thus termed stochastical/statistical estimation averaging, more generally, we call it estimation averaging or filtering. The estimation averaging may be performed on the grid-scale of the engineering model, o o Y(x)s = j H(x - x')p H(x - x')s X(x') dx', (1.5) —oo where Ys is the estimated X in the whole domain and Hs is the filter epitomizing the 5Here the closure model refers to the one that does not require an explicit description of the dynamics on scale smaller than Ap [Beckie, Aldama and Wood, 1994]. In the following discussion, the 'measure' implies what is meant by Cushman [1986]. 6 T h e variable state in the unsampled region is somehow an average of the variable state measured in the sampled region. In the case of unconditional simulation, it still bears an essence of averaging because the inferred model structure is established in an averaging mean sense. Chapter 1. INTRODUCTION 7 estimation procedure. If on a scale smaller than a grid, a block averaging is then applied, oo Y{x)B = j H(x- x')p H{x - x')s H{x - x')B X{x') dx', (1.6) —oo where YB denotes the averaged block estimates of X, and HB is block averaging filter. This process belongs to the so-called upscaling, which shall be discussed in next section. Now we see that the estimation averaging always invokes measurement averaging, thus, its accuracy depends upon the accuracy of the measurement averaging. This fact seems frequently be ignored in the literature. For example, on the conditional simulation work of Delhomme [1979], two assumptions have been made: 1) the simulated values of conductivity have the same correlation structure as the 'true' conductivity field, and 2) the simulated conductivity at the measurement(e.g., pumping test) location are consistent with the 'observed' values. It is clear that the so-called 'true' or 'observed' conductivity is inverted from the drawdown recorded at pumping and observation well [Theis, 1935; Cooper and Jacob, 1946], representing the spatial/spatial-time averaged aquifer property through the associated measurement filters. Thus, the inverted conductivity is neither 'true' nor 'observed', instead, an averaged value. Consequently the accuracy of the estimation averaging is limited by the accuracy of measurement averaging. Now it becomes clear that the engineering world only possesses finite scales, not only because the measurement process can only sample scales of limited ranges, but the dy-namics become invalid beyond certain scale — the resolving scale. For example, Darcy's law is not valid on molecular scales and remains unproved in large-scale heterogeneities. The effect of unresolved small-scale dynamics on resolved large-scale dynamics can be lumped into phenomenological parameters on the resolving scale. For example, molecular process — • fluid continuum — • porous media continuum energy conservation p = ^ jf- (Newtonian fluid) K = ^ k Chapter 1. INTRODUCTION 8 where // is viscosity [ M L - 1 T _ 1 ] , c for tangential stress[ML - 1T - 2] and jj"- represents veloc-ity gradient. The phenomenological parameter fi lumps the molecular moments c; and while the conductivity K integrates the effects of the fluid ^ and the pore k. Indeed, it has been proved that Darcy's law can be directly derived by averaging the Navier-Stokes equations Keller [1980]. Beven [1989] recognized the effects of the interaction between the subgrid-scale (un-resolved) and the grid-scale(resolved) on the resolved dynamics, in the evaluation of the physically based models, and he concluded that the present physically based models essentially still belong to traditional lumped conceptual models. The subgrid-scale ~ grid-scale interactions in groundwater flow models have been recently, for the first time, evaluated by Beckie, Aldama and Wood [1994]. These authors adopted some averaging concepts from turbulence continuum and statistic physics literature and employed a spa-tial filtering approach(references therein). They found that if the resolved and subgrid scales of the variability of the porous medium are separated by a spectral gap, Darcy's law manifests an accurate resolved scale model structure. The so-called spectral gap A / equal to c5 — 7 should meet the condition, A / > 2 7, (1.7) where in frequency space, [—7,7] represents the resolved scale range(or wavelength equal to / _ 1 ) , and [—fi, — 8] and [8, u] denote the subgrid scale ranges on two sides of zero frequency. The preceding paragraphs complete the discussion of Fig. 1.1. To further our under-standing, we graphically show in Fig. 1.2 the scale and engineering averaging relationships between real world, measurement, model and decision analysis. The two circles represent measurement filtering, and estimation filtering and model upscaling. Note that these graphs demonstrate scale versus correlation range (A ~ a), and they differ from the REV Chapter 1. INTRODUCTION 9 A Real World JL C Measurements soft-data-technique ^remote sensing] tracer test -y— core test well test Up E Engineering Model Figure 1.2: The A(scale) ~ a(correlation range) of continuously filtered real world(A), measurements(C), model(E), decision analysis(F), and filtering processes(B, D). Chapter 1. INTRODUCTION 10 illustrations of Bear [1972] (Fig.1.3.2), Cushman [1984] (Fig.l), and Baveye and Sposito [1984] (Fig.l). What the authors show is averaged parameter versus the scale of aver-aging region or volume. Here, Fig. 1.1 A manifests the conjecture that correlation range continuously increases with real world scales, but decreasing f j ^ j . This can be best explained by a hypothetical filtering [Vanmarcke, 1983]. Consider an infinite continuous signal, we filter this signal by a bell-shaped filter. We then can see an increasing correla-tion in the filtered signal from a filtering with larger filter width and this increase slows down with the further increase of filter width. Then graph A is a plot of the filter width versus the correlation length of the filtered signal. Fig. 1.1C shows the measurement filter width Ap versus the correlation range of the sampled signal. Here we have assigned a representative scale to each measurement, accordingly, we observe the same scale relationship as that in graph A except here both Ap and ap are in discrete form. Note that the 'soft-data-technique' may be able to achieve 'soft' continuous relationships with 'experience filters'. Also note that ap doesn't start from the origin point because the measurements can not sample the scale smaller than Ap, thus the correlation is uncertain below some cut-off value. This differs from the correlation range of a conventional variogram of point data support. The model resolving scale can be represented by its grid block scale, A#. Shown in Fig. 1.1E, the A^ and aB of a model are relatively more uniform than those of the mea-surements, the decision variables and the real world. The scale of the decision variable, shown in Fig. 1.1F, Ap, varies for different decision variables and for different engineering applications. The realistic Xry and aD, being concerned by the decision maker, should be consistent with the Ap and ap of the measurement and the \B and aB of the model. In other words, the measurement should be designed to have the Ap and ap consistent with \D and aD, so does the model. Chapter 1. INTRODUCTION 11 Meanwhile, model dynamics should be scaled up. In water resources context, dy-namics upscaling has been pioneered by several authors, Hassanizadeh and Gary [1979a, b]; Baveye and Sposito [1984]; Cushman [1984, 1986]; and Beckie, Aldama and Wood [1994]. The basic idea employs a filtering approach, by applying appropriate filters to filter the small scale mass, momentum, energy and entropy balance equations, upscaling to a large macroscopic scale amenable for field measurement and model representation. These filters can be any or all of the measurement filter Hp, estimation filter Hs, and block scaleup filter HB. Dynamics upscaling is beyond the scope of the thesis, thus it shall not be discussed further. Parameter upscaling is discussed in the following. 1.1.3 Parameter upscaling: A filtering interpretation Filtering from a smaller scale to a larger scale, through the filters Hp, Hs, HB, is ex-actly the essence of upscaling. Parameter upscaling is an engineering compromise. The measurement filter Hp is utilized as an instrumental technique scaling-up and scanning-out unmeasurable real-world small-scale variability. On the other hand, the estimation filter Hs, upscaling to the whole domain, is a practical concession to cost and techni-cal constraints. Furthermore, if the block-scale of the engineering model is greater than the data support of measurement and estimation filtering, the block scaleup filter HB is employed to provide data of the model-input format. Consequently, upscaling concep-tually unifies effective parameter, equivalent parameter and measurement process, since all are scale-dependent concepts and related to some basic small-scale-unit, e.g., core permeability. In Appendix A, centering on upscaling we review the work on effective and equivalent parameter, in particular on conductivity, as a limiting case and a general case of upscaling, respectively. Chapter 1. INTRODUCTION 12 1.1.3.1 Parameter measurement: a special case of upscaling The measured parameter can be considered as a special case of upscaling, in the sense that the measured parameter is a result of the upscaling through technical measurement. In this section we specifically review the work on pumping tests in relating measured conductivity to the small-scale conductivities. We denote Ym the logarithm of the mea-sured conductivity, and Y the logarithm of the small-scale conductivity. In next section we discuss the work on slug test — the theme of the thesis. In the work on pumping test in confined aquifer, Butler [1991] found that Ym is dependent on the spatial location of observation well employed in the analysis methods of Theis [1935] and Cooper and Jacob [1946].7 In the analysis of angular dependence, he found that the variability of Ym increases with try-(variance) and the distance from pumping well, so is the radial dependence, but the radial dependence is much greater than the angular dependence. He stated that the dependence should be weakened with decreasing correlation of Y. Using a semi-analytical method, Butler and Liu [1991] worked on pumping test sited on the center of an infinite linear strip of one material embedded in a matrix of differing hydraulic properties. The drawdown contours are shown for several cases: Kst"P/Kmatrix=0.1, 100, 10000(see their Fig. 5, 7, 8)8. Sensitivity analysis is employed to investigate the influence time of the strip on drawdown, depending on strip width and the Kstr,pI'KmatTtx ratio. Various settings of pumping and observation wells relative to the strip location are discussed. It is found that the strip can only affect the changes in drawdown for a finite time, meaning that after certain time further changes will be independent of K3trxp, although the total drawdown will always be dependent on Kstrtp. 7 T h e spatial dependence includes angular and radial dependence which refer the direction and the distance of observation well to pumping well, respectively. 8 K denotes conductivity. Chapter 1. INTRODUCTION 13 This is similar to the sensitivities of zones lying radially outward from the pumping well and observation wells of Butler and McElwee [1990]. Butler [1990] found that the log-log [Theis, 1935] and semilog [Cooper and Jacob, 1946] methods yield dissimilar estimates due to their emphasis on properties in different portions of the aquifer. The difference is proportional to Cy, and inversely proportional to the outward distance from pumping well if observation well record is used. Butler [1990] discussed two interpretations: (i) Using pumping well drawdown: the near-well conductivity strongly impacts the K from the log-log method^"), whereas the K from semilog method (Ksl), using large time portion, can not "feel" the effect of the near-well conductivity, and it represents the property of the aquifer at the front of depression zone, (ii) Using observation well drawdown: the effect of the near-well conductivity on Kn lessens with increasing outward distance from pumping well. Furthermore, Butler [1990] stated that the close agreement between the slug test and large-time Ksl at several observation well should be considered as strong support for representing the aquifer as a uniform at the slug-test or larger scale. One crucial point Butler [1990] did not explicitly elucidate is: the dissimilar estimates of Ku or Ksl, at the pumping well and the observation well, result from the differing flow patterns around the pumping and observation well, i.e., the filtering power and volume around the pumping and observation well are different [Oliver, 1993]. Butler and McElwee [1990] utilized the sensitivity analyses of observation well draw-down to T and S, denoted Uj and Us- By the analysis of sensitivity in the zonated aquifers simulated numerically, they found that, only during a finite interval of time, UT is sensitive to the near-pumping well material, being consistent with the analytical resu\ts[Barker and Herbert, 1982; Butler, 1988]. This is also true for Us- They concluded that drawdown at an observation well placed at a distance from pumping well is rela-tively insensitive to material lying between the pumping well and the observation well, as Chapter 1. INTRODUCTION 14 only the early-time responses are dependent on the properties of the interwell material. Consequently, observation wells are superior to more distant wells for the purposes of characterizing property variability in a radial-flow field. Their findings can be examined by the filters of the pumping and observation wells [Oliver, 1993]. In a uniform aquifer, McElwee and Yukler [1978] analytically derived the sensitivities: Ur = ~ + Sfi exp (- ^ ) (1.8) T 4irT2 -Q 4TTT5 Us = TS-C exp (^) (1.9) 0 0 -u s = Q f e j €—du (1.10) 4TTT r*S/ATt One interesting point can be drawn from the formulas is that systems with lower T, will have higher sensitivities and thus will, in general, be more amenable to inversion analyses, given a sufficient time of pumpage has elapsed [McElwee, 1987]. Oliver [1990] characterized the relationship between drawdown and permeability as a Frechet kernel function of the radial flow in pumping tests, Si = / . G , ( x ) ( l - ^ y ) 4 ; KD(x) = (1.11) where Gi(x) is the kernel and (1 — K^x^) is the presumed permeability distribution which is radially symmetric but nonuniform, and K is the average permeability. He solved for the kernel by a perturbation method. This kernel function was used by Oliver [1992] to derive the permeability distribution of eq.(l.ll). Oliver [1993] generalized the (1 — K^x^) permeability pattern to nonradially symmetric porous media and analytically solved the associated Frechet kernel by a perturbation method. Desbarats [1994] defines the effective conductivity measured from 3D flow pumping test analysis as a weighted "power average", Chapter 1. INTRODUCTION 15 where V is drainage volume. The | ^ j r weight was first seen in the work of Cardwell and Parsons [1945] in the determination of the upper and lower bounds of Kev- Desbarats [1992b] also employs the f^jp weight in the definition of Key in 2D radial flow. The power exponent w in eq.(1.12) is obtained by fitting Kev to numerical simulations. 1.2 T h e measurement filter: slug test In Section 1.2.1 we first review the slug test literature, then we discuss the slug test from a spatial filtering perspective in Section 1.2.2. 1.2.1 Review of slug test literature Hvoslev [1951] may have been the first to point out that a slug test can be used to estimate in situ hydraulic conductivity. Hvoslev [1951] develops models for a variety of well, piezometer, and standpipe geometries with full or partial penetration in an infinite, semi-infinite below, or confined homogeneous, anisotropic aquifer. The assumption of no aquifer storage in the models of Hvoslev [1951] was evaluated by Chirlin [1989]. Demir and Narasimhan [1994] developed an improved interpretation strategy of the slug tests of Hvoslev [1951]. Cooper et al. [1967] presents the mathematical model of the slug test of a finite wellbore radius in infinite confined homogeneous media. By the analogy to the heat conduction problem of Carslaw and Jaeger [1959], Cooper et al. [1967] give the analytical solution of hydraulic head at wellbore and in aquifer. Papadopulos et. al [1973] extended the type curves of Cooper et al. [1967] to those with Sa = 10 _ n, n = 6, . . . , 10 [L - 1 ] . 9 9 For simplicity, the dimension unit of specific storage is hereafter omitted. In other words, we consider an aquifer of unit thickness. Chapter 1. INTRODUCTION 16 Sageev [1986] discussed the influence of wellbore skin and dimensionless wellbore storage10 on the early-time, intermediate time and late time type curves. Karasaki et. al [1988] presented analytical solutions of slug tests in bounded porous media, fractured rocks, with linear and radial flow. Their results show that slug tests suffer problems of nonuniqueness to a greater extent than other well tests. Guyonnet et. al [1993] empirically determined the distance traveled by a given head perturbation H/HQ through regression analysis for the slug test in unbounded homoge-neous media. For H/Ha = 1%, where r(t) the traveled by the H/H0 = 1% perturbation at time t, and r m a x is the maximum traveled distance. In the analysis of conductivity measured by a slug test in confined heterogeneous conductivity field with square blocks on wellbore scale, Harvey [1992] employ the power-average definition Kev of eq.(A.II.152) and determine the exponent w by fitting Kev to the results of Monto Carlo simulations. The integral of eq.(A.11.152) is evaluated by summing for all conductivity block values within a cylinder of radius ropt with power w°pt where ropt and wopt are the values that provide best fit of Kev to numerical results. Harvey [1992] states that (i) ropt linearly varies with (ii) wopt decreases with xJjbS , where b is aquifer thickness, rw wellbore radius, Ss specific storage, and \h horizontal correlation. Beckie and Wang [1994] developed a spatial filtering approach to relate the conduc-tivity measured by a slug test to the core-scale conductivities(on wellbore scale). This 10Wellbore storage is 2*r3"'i 5 , where S is storativity and re and rw are well casing radius and well screen radius. (1.13) (1.14) Chapter 1. INTRODUCTION 17 shall be discussed in the next section. 1.2.2 The slug test filter Unlike any methods reviewed above in determining the effective or equivalent parameter or the parameter measured from aquifer tests, the spatial filtering method explicitly links the core-scale conductivity to the measured conductivity through a spatial filter, which quantifies the contributions of core-scale conductivity at different locations to the measured conductivity, where Yc(x), Ym(x) are the logarithm conductivities of the core-scale field and measured field, and G(x) is the spatial filter. The filtering representation and filter identification are the subjects of the following chapters. In order to obtain a stable slug test filter,11 the spatial filtering approach requires a great number of measurements, and the more the measurements, the better the stability of the filter. In the thesis, the number of slug tests is at the magnitude of 103 ~ 105. Note that the number of measurements is far more than the number of Monto Carlo runs conducted by various workers [Harvey, 1992; Desbarats, 1994; Butler, 1990], who intend to relate the statistics of Ym to that of Yc. Thus, the numerical spatial filtering approach may provide results of superior stability. On the other hand, the noise in the measurement process can be evaluated and removed through an optimal Wiener filtering technique, whereas it would pose difficulties for the Monto Carlo method. Furthermore, the spatial filtering method is versatile to account for conductivity field of high vari-ance, multimodal distribution(e.g., bimodal sand-shale) and nonstationarity and nonlin-ear transient flow, whereas these conditions are hard to be accommodated by analytical 1 1 Here 'stable' means that the adding of new measurements will not improve or degrade the identified slug test filter. (1.15) Chapter 1. INTRODUCTION 18 methods(e.g., first-order perturbation, spectral method) without simplifications. The en-gineering application and theoretical significance of the spatial filtering approach shall be summarized in Chapter 7. 1.3 Thesis overview The thesis subject is the spatial filtering from core-scale conductivity to the conductivity measured by a slug test and the effect of specific storage on the conductivity filtering through numeric experiments. In Chapter 2, we present numerical methods to simulate slug test in heterogeneous and homogeneous media. The conceptualization of slug test in a system and filtering framework is given in Chapter 3, and in which system filtering properties and the inverse estimation of system outputs are discussed in detail. In Chap-ter 4, a spectral analysis through Wiener filtering is developed to identify spatial filter in heterogeneous media. As a result of computation load constraint, the specific storage value examined in the heterogeneous conductivity filtering is limited to Ss = 10 - n , n = 1, . . . , 5. The effects of other S3 with values of n = 6, . . . , 10 on the conductivity filtering are investigated in homogeneous media where the spectral analysis method becomes invalid, a numerical perturbation method is developed, which is the topic of Chapter 5. In Chap-ter 5, empirically deterministic expressions of influence radius and equivalent filter width versus specific storage are presented. Consequently, concerning the aquifer considered in Chapter 4 and 5, the two chapters contain conceptually different but complementary ap-proaches. In Chapter 6, the conductivity filtering in heterogeneous media with different specific storages and various spatial variabilities are explored. In Chapter 7, the major results are summarized in conclusions and future work is recommended, also the possible engineering application and theoretical significance of the spatial filtering approach are discussed. Chapter 2 S I M U L A T I O N O F S L U G T E S T In this chapter, we introduce the boundary value problem associated with the slug test and its numerical solution. For the research in subsequent chapters, we are required to simulate the slug tests in the following domains: a) Two dimensional, heterogeneous conductivity field; b) Two dimensional, homogeneous conductivity field with a small region of inhomogeneity, or a so-called perturbed homogeneous field. c) One dimensional, homogeneous conductivity field. We utilize three different finite difference models to simulate slug tests. Considering the constraints in the simulation, i.e., 1) the statistical requirement of random conductivity, 2) the hydraulic gradient distribution excited by the slug test, and 3) the efficiency, we discretize the domains with three appropriate type of grids, shown graphically in Fig. 2.3. Details of these discretizations will be discussed in the corresponding sections. This chapter proceeds as follows: First, we present the slug test set-up and its an-alytical solution. Second, we introduce the three FDM models, the approximation and discretizations techniques. For illustration convenience, we first discuss the 2D hetero-geneous case, then the 2D perturbed homogeneous case, and lastly the ID homogeneous simulation. 19 Chapter 2. SIMULATION OF SLUG TEST 20 © ID homogeneous field radial symmetric coordinates Figure 2.3: Domain discretization of slug tests. Graph A is a rectangular discretization in heterogeneous conductivity field; B is a radial discretization in homogeneous field with a small region of inhomogeneity(only half shown); C is a radially symmetric discretization in homogeneous field. Chapter 2. SIMULATION OF SLUG TEST 21 2.1 T h e boundary value problem of a slug test Fig. 2.4 shows the schematic diagram of a slug test of a finite wellbore diameter in a infinite confined aquifer. The flow, initial and boundary conditions will be described by a mathematical model after the next paragraph. / / / / / / / / / ? / / / / / / / / / Figure 2.4: A slug test in a confined aquifer. To facilitate the numerical simulation of above problem in heterogeneous media, we simplify the problem by assuming that the vertical flow is negligible, by considering that the wellbore is screened throughout the whole thickness of the aquifer, and that the aquifer stretches to infinity. It is then not unreasonable to assume the horizontal hydraulic gradient to be much greater than the vertical hydraulic gradient, i.e., | £ >• | j . Consequently, the assumption of negligible vertical flow has the effect to reduce the flow dimension from 3D to 2D. Other practical issues are excluded from this thesis, e.g., the effects of initial excess head and wellbore damage, the double straight line and inertial effects, and the effects of entrapped air and oscillatory responses etc.. Chapter 2. SIMULATION OF SLUG TEST 22 The simplified model is represented mathematically, V-(K-Vh) = Ss^, r > rw, (2.16) K—rdOdr = *r2c — , r = rw, (2.17) o h(rw,t) = Hit), (2.18) H = H0, t = t0, r = rw, (2.19) h(r,0) = 0, r > rw, (2.20) where, K: conductivity. Sa: specific storage. Ho: initial excess head. H(t): water level at wellbore at time t. h(r,t): head in testing aquifer at distance r,time t. rc: well casing inner radius. rw: well screen inner radius. Note that, by the principle of superposition, the initial conditions can be generalized to any case. Cooper et al. [1967] presented an analytical solution of a slug test in homogeneous medium. This solution was taken from heat conduction literature [Carslaw and Jaeger, 1959]. Hydraulic head at wellbore is solved as *(,) - i i l ^ / V * * - * (2.21) TT2 Jo U • Au Au = [u JQ(u) - 2 a J1(u)]2 + [uY0(u)-2 aYiiu)]2, a = — ss, a = — — , rc rc Chapter 2. SIMULATION OF SLUG TEST 23 where S3 and K represent specific storage and conductivity, respectively; and ct and B are of the same dimension [L - 1]. The so-called dimensionless wellbore storage [Carslaw and 2 Jaeger, 1959; Sageev, 1986; Karaski et ai, 1988; Guyonnet et ai, 1993] is 2 * r Z c s = jsJH' where S is storativity and b is aquifer thickness. In the thesis we assume rc = rw, and use a different terminology. We call the a = Ss at wellbore the wellbore specific storage, or for short, wellbore storage. Because the wellbore is a hollow cylinder only containing water, thus wellbore storage equals to 1 [L - 1 ] . 1 For slug test in the simplified 2D heterogeneous aquifer, the solutions can be obtained numerically. In the next section, we will introduce a FDM model in such an aquifer and discuss in detail the numerical approximation errors and techniques. 2.2 2D F D M in rectangular coordinates For the work that follows, we must simulate slug tests in heterogeneous conductivity fields. The conductivity fields are generated by geostatistical simulation methods which utilize regular grids. In numerical models, the fields are most often discretized in blocks regularly spaced in cartesian coordinates. For rectangular grids with constant mesh spacing, finite difference and finite element methods have no advantage over each other. In the thesis, we select the finite difference method because of its appealing simplicity in discretization and the amenability for an efficient symmetric matrix. The derivation of the FDM equation is depicted in Appendix B.I. Here only the aspects affecting the accuracy of the FDM method are discussed. 1 For simplicity, the dimension unit of specific storage is hereafter omitted. In other words, we consider an aquifer of unit thickness. Chapter 2. SIMULATION OF SLUG TEST 24 2.2.1 Numerical approximation and discretization errors The discrete output of numeric models are utilized to approximate the continuous re-sponse of the physical systems. It is well known that that we can use the sampling the-orem to reproduce a continuous function through discrete sampling [Oppenheim, 1983], hence it is possible to use the sampling theorem as a guide for the time and spatial dis-cretization of numeric models. However, this is very hard to achieve when the physical response is governed by a set of intractable partial differential equations predicated on the boundary and initial conditions. Thus, here we follow a 'try-evaluate-try' methodol-ogy which rests on the proper physical interpretation of model outputs. We minimize or exclude the time and spatial approximation errors until there is a best fit between the model outputs and the expected outputs possibly obtained with the accurate analytical solutions of some idealized systems, which in this case is the slug test in homogeneous media. Even though the slug test is one of the simplest hydrogeologic tests, its numeric simulation is far from a trivial task. The complexity arises from the extremely dynamical response and the Dirac initial conditions of hydraulic head across the wellbore and the interactive role of conductivity and specific storage in the recovery. These aspects are discussed in the next 4 subsections and a summary of numerical errors is given at the end. 2.2.1.A Wellbore approximation In the 2D FDM model, the wellbore is represented with one grid block to facilitate computer programming(see Fig. 2.3 A). We numerically derived in Appendix B.I that the wellbore grid is equivalent to a real well by assigning wellbore storage 5^e" = 1 and wellbore conductivity Kwelt several magnitudes higher than that in the aquifer. Chapter 2. SIMULATION OF SLUG TEST 25 S,""1 = 1 K « U / K ^ u l f c r > 1 ( ) S \ ) Wellbore 0 * w2< / A OZ. 2L, i i r. N 1 / r K apparent screen radius Ah A2P2 A h " 2 " 2 A h A ' P l Ah (A) u u u i i / \ i 1 U I I Ax Ax s Ax ^ Ax (B) W j P l Figure 2.5: A: Discretizing wellbore into a number of fine blocks. B: A round well is approximated by a square block. This treatment of the wellbore is particularly important for the simulation of slug tests. In this case, a general discretization formula can be applied at the wellbore with the only need of assigning S™el1 and KwM to the wellbore grid. More importantly, it is feasible to assign any number of blocks as wellbore so as to readily achieve a finer wellbore discretization that is shown in Fig. 2.5 A, or, investigate the spatial discretization errors, which will be discussed in the next section. Approximating a round wellbore with radius ra by a square grid block with dimension Ax related to ra by 7T r 2 = Ax 2 , ra = 0.564 Ax, (2.22) will introduce numerical error even though it implements the mass balance equation (2.17) exactly. Shown in Fig. 2.5 B, it overestimates the hydraulic gradient | £ W i P 2 while un-derestimating ^ W l P l . This shape effect becomes less significant for a finer discretization Chapter 2. SIMULATION OF SLUG TEST 26 proximal to wellbore. 2.2.1.B Spatial discretization Hi Ho Hk head distribution at time to head response at wellbore aquifer initial head I ° rw r ° t0 t (A) (B) Figure 2.6: Head distribution at initial time and head response at wellbore. Figure 2.6 A shows head distribution at initial time. The Dirac-function-like feature of the head distribution across the inner boundary requires the space discretization proximal to wellbore to be infinitely small at the very early time. As time progresses, this requirement eases. The resultant numeric error can be reduced by refining the spatial discretization, simply by assigning a number of blocks at wellbore and employing a fine discretization in the vicinity of wellbore. Fig. 2.7 shows the simulated response curves of two discretization schemes denoted as a and b. Curve a represents the discretization of evenly spaced grids of a wellbore diameter(hereafter denoted as Axw), while curve b stands for the grids of 1/13 of a Axw. We see that curve b matches the analytical solution very closely. This confirms that the misfit of curve a is a result of coarse spatial discretization. The type of spatial discretization error shown by curve a can also be observed in the work of Demir Chapter 2. SIMULATION OF SLUG TEST 27 Dimensionless time Figure 2.7: Wellbore response comparison: numerical solution against analytical solution, a: space step equal a wellbore diameter; b: space step equal ^ of a wellbore diameter. All other parameters of a and b are identical. and Narasimhan [1994] who utilize an integrated finite difference method to simulate slug test and treat the wellbore as a capacitor (see their Fig. 5). The Sa of Fig. 2.7 is 10 - 1. With other Ss values, Fig. 2.8 displays the close match between the analytical solutions and the numerical results of refined meshes. The evenly spaced grids are of 1/5, 1/4, 1/3, 1/2 of a Axw, for Sa equals to IO - 2 , IO - 3 , 10~4, IO - 5 , respectively. For Ss < 10~5, the numerical errors of the grids of a Ax*" become negligible. This means the spatial discretization error decreases with decreasing Ss. It is because a small Sa invokes a large influence zone, hence small hydraulic gradient, and mitigates Chapter 2. SIMULATION OF SLUG TEST 28 Dimensionless time Figure 2.8: Wellbore response comparison: numerical solution against analytical solu-tion. The solid curves are analytical solutions; The dashed curves are simulated results. Specific storage ranges from 10"4(upper right) to 10_ 1 (lower left). Their space steps are | , | , of a wellbore diameter respectively. spatial discretization. However, such grid refinements increase the number of degrees of freedom significantly. For example, with Ss = 10 _ 1, and considering the simulation of a domain of 64 x 64 Axw, the number of freedom after refining the grid to 1/13 of a Axw will be 64 x 64 x 13 x 13 = 692224 blocks. For smaller 5S, although using a large grid block of a Axw may be feasible, the drastic increase of flow domain requires overwhelming grid size(the number of grids). Computation limitations suggest it is better to use the spatial discretization of one grid block for the wellbore. The resultant numerical errors can be mostly offset by employing Chapter 2. SIMULATION OF SLUG TEST 29 an effective radius which is introduced next. We define the effective radius as that which yields the best match between a numer-ical simulation in a homogeneous aquifer and the analytical results [Cooper et al., 1967; Beckie and Wang, 1994]. Beckie and Wang [1994] found r e = 0.59Axw for S8 < 0.15, while Harvey [1992] set re = 0.427Ax™ and Peaceman [1978] used re = 0.19Axu'. The dif-ferences between these workers arise from the methods used to approximate the wellbore. Logrithmic specific storage Figure 2.9: Effective screen radii in terms of wellbore diameter versus logarithmic specific storages. Figure 2.9 illustrates r e versus Ss. We see that all effective radii are greater than the apparent screen radius r„ = 0.564Ax*".2 Using re > rw reflects the attempt to increase flux. On the other hand, it implies that in general the numerical errors inhibit flux to the wellbore. This is true because the approximation of | £ by ^ underestimates the spatial gradient, particularly near the wellbore. We also see that, for Ss smaller than 2Axw = Ax, if the wellbore contains only one grid block. Chapter 2. SIMULATION OF SLUG TEST 30 10~6, the effective radii approaches the apparent screen radius. This is consistent with our observation that the spatial discretization errors are negligible for Sa < 10 - 5. The numerical spatial discretization error varies with time. The time discretization will be discussed in the next section. Here we only show that if we evaluate re at each time step, we end up with a time-dependent effective radius re(t). i — i — i — r re(t)[Ss=l.e-l] -O-re(t)[Ss=l.e-2J re(t)[Ss=l.e-3] re(t)[Ss=l.e-4] • re(t)[Ss=l.e-51 • 9 13 17 21 25 29 33 37 41 45 49 53 Time steps Figure 2.10: Effective screen radii in terms of wellbore diameter (Ax10) are a function of time. The smaller Sa, the flatter the re(t), the faster re(t) reach 0.564Axto. (The last time step is, t55 = ^ = 1 M 2 ^ ! 1 ! ) . Figure 2.10 demonstrates the time dependence of the estimated re(t), for 10 - 5 < Sa < 10_ 1. We see that re(t) curves approach the apparent well screen radius 0.564Axw at later time steps, when the hydraulic gradient becomes flat enough to assuage discretization error. We also see that the smaller specific storages, the faster re(t) falls off, and the closer re(t) stays to the apparent well screen radius 0.564Ax™. For Sa < 10 - 5, the re(t) curves overlap the 0.564Axtw line. The reason for this is again the effect of Ss on the Chapter 2. SIMULATION OF SLUG TEST 31 hydraulic gradient distribution. 2.2.1.C Time discretization We have shown, in Figure 2.6 B, the wellbore response with respect to time. Starting from time to, the response is continuous, this fact eases the time discretization. The static hydraulic heads are almost recovered as 8 = approaches 100 [Papadopulos, 1973] and the recovery becomes independent of Ss-Considering that the recovery becomes insignificant at later time for Ss > IO - 5 , we simulate our slug tests in heterogeneous media until 8 = 10. The total real time tmax is calculated by w (2.23) "max — Ty 5 where KQ is the geometric mean of the heterogeneous core-scale conductivity. We dis-cretize tmax to steps geometrically increasing by a constant factor tf. Thus the first time step t\ should be 4 i - i - (tfr ( } where nt is number of time steps. If 8 and t j are specified, the t\ depends only on nt. The greater the specific storage, the faster response, thus the more time steps are required. For a given Ss, we simply adjust nt until the time discretization error dissipates. In the thesis we choose nt — 55 for tf = 1.15 or nt = 100 for tf = 1.1, and we find that time discretization error becomes negligible for IO - 1 0 < Sa < 10_ 1. 2.2.1.D The hypothetical remote boundary Mathematically, the model solution requires a remote boundary condition at an infinite distance from the wellbore. The effect of a finite domain be minimized if we place the Chapter 2. SIMULATION OF SLUG TEST 32 boundary sufficiently far away from the wellbore. However, this results in an excessive computational burden because the grids are evenly spaced and of finite length. To get around this, an approximation technique is employed to divide the whole domain into a heterogeneous inner zone with blocks equally spaced and a homogeneous outer zone with blocks of length increasing by a geometric factor each step away from the inner zone [Harvey, 1992]. The homogeneous conductivity is set to the geometric average of the inner zone. The diagram was shown in Fig. 2.3 A. The concept of influence zone is frequently seen in the well testing literature, yet no universal definition can be found. It is due to the 'dynamic' essence of what is meant by 'influence'. For slug tests, for example, it is called "range" of the slug test [Barker and Black, 1983]; "radius of investigation" [Sageev, 1986]; "the maximum distance traveled by a perturbation jjr = 1%, 5%, 10%" [Guyonnet et ai, 1993]. A more rigors definition will be given in Chapter 4, here we temporarily adopt the -jg^ = 1% definition of Guyonnet et al. [1993]. The influence zone at large enough time(e.g., Q = 100) is determined by Ss and is independent of K. Guyonnet [1993] experimentally determined that area of influence zone is proportional to specific storage, r/rw oc [xrl/ftwrl, b S^)]0-495. Our initial results show that, for I O - 1 0 < Sa < 10 - 1, 20 blocks is enough to exclude the boundary effect if the geometric factor is set to 2. The number of blocks at inner zone is obtained by enlarging the inner zone until wellbore response is not affected by the presence of the outer homogeneous zone. It is found that 64 by 64 blocks is sufficient for 10~3 < Ss < 10 - 1. More blocks are needed for smaller Ss and computation becomes expensive. In the preceding discussions, we only used the wellbore response curve to investigate the numerical errors. This is because wellbore response is most dynamical so that its simulation possesses the maximum error. We are assured that the simulation error within flow domain should be trivial if the simulated wellbore response curve is accurate enough Chapter 2. SIMULATION OF SL UG TEST 33 to have matched the analytical solution. In summary, the sources of numerical errors include: (a) The temporal and spatial discretization by approximating | £ | | with ^f, Using a /('-dependent discretiza-tion, the temporal approximation can be reduced to a negligible level. But the spatial discretization error can not be easily eliminated due to computation constraint, (b) The boundary approximation. The wellbore shape effect is unavoidable but not significant, while the effect of the remote boundary can be minimized. Errors of type (a) consistently cause an underestimation of the hydraulic gradient, which can be offset by applying a effective radius greater than the apparent radius. Errors of type (b), for approximation of a round wellbore, cause under- or overestimation of hydraulic gradient, while for the remote boundary, it may positively contribute the wellbore flux if not placed far away. 2.3 2D F D M in radial coordinates Simulation of a slug test in homogeneous confined aquifer with a region of inhomogene-ity(Fig. 2.3 B) can be conducted by a 2D FDM model in radial coordinates. The radially distributed mesh conforms with the flow pattern invoked by a slug test and thus is much more resistant to space discretization errors, and is more efficient in reducing the number of degrees of freedoms. The mathematical model is exactly the same as that of the slug test in heterogeneous aquifer introduced in the previous section. We only need transform the heterogeneous model to radial coordinates and perform the same finite differencing. We rewrite the flow equation (2.16) as The detailed FDM equations are derived in Appendix 2.3. In the FDM model, besides the boundaries of the 2D heterogeneous model, one more Chapter 2. SIMULATION OF SLUG TEST 34 boundary of 0 is needed. We set the 9 boundary to be opposite to the perturbed inhomo-geneity and treat it as a nonflow boundary. Note that there is only one inhomogeneity and it is placed along a centerline of the grid. The nonflow boundary is physically valid because the symmetric flow pattern invoked by the inhomogeneity warrants no crossing flux. Note that we can not simply simulate only some slices of the radial domain and incorporate the 0 boundary effects by means of superposition like that of pumping wells. Because the boundary effects on the inhomogeneity affecting the flux to the slug well differs from the boundary effects on the pumping or of any other source terms. Thus, we need simulate the whole flow domain, which was shown in Fig. 2.3 B. Note that only half of the domain is illustrated. The spatial discretization, unlike the constant grid spacing in the 2D heterogeneous model, can be easily refined until sufficient accuracy is achieved. The | £ discretization is obtained empirically as a function of Sa. The | | discretization of 50 equal A0 is found to sufficient. The time discretization formula is identical to that of the 2D heterogeneous model, but for better accuracy 100 time steps and a factor of 1.1 are used. Apparently, the wellbore shape effect does not exist in the radial model. Fig. 2.11 demonstrates the excellent agreement between the numerically simulated wellbore response and the analytical results. The minor graphical misfit is due to the inadequate data points in the original paper [Cooper et al., 1967], 2.4 I D F D M in radial coordinates We use a ID FDM homogeneous model to reproduce the analytical solution of Cooper et al. [1967] and Papadopulos et al. [1973]. The ID FDM model is very important to efficiently reproduce the wellbore response curve needed in the inverse estimation of the measured conductivity discussed in Chapter 3. Even though the analytical solution Chapter 2. SIMULATION OF SLUG TEST 35 Dimensionless time Figure 2.11: Simulated wellbore response(dashed lines) versus analytical solution(solid lines). The graph only shows the curves of Sa ranging from 10 - 5 (upper right) to IO - 1 (lower left). For Ss < 10"5, the fitting improves. can be evaluated by the numerical inversion of Laplace transform with the algorithm of Stehfest [1970], which is found to be of great use in hydrologic applications[Moenc/i and Ogata, 1984; Butler and Liu, 1991; Guyonnet et al., 1993; Harvey, 1992], the efficiency of this algorithm is still not adequate for the problem at hand. For example, considering 256x256 slug tests are recorded with each comprising 55 time points, and assume the inverse fitting procedure(discussed in Chapter 3) converges by 10 searches, then the number of calls to the Stehfesfs [1970] algorithm is 3.6 xlO 7, this is overwhelming! The example is just a moderate case of those discussed in Chapter 3. Chapter 2. SIMULATION OF SLUG TEST 36 However, the homogeneous model can be easily solved by the ID FDM. The straight-forward boundary conditions can be readily implemented. The spatial discretization is identical to the | £ discretization of the 2D perturbed homogeneous model, and the time discretization is also identical. We rewrite the flow equation (2.16) in a one dimensional form, 1 d dh Sa dh rcTr{rTr) = K Tf (2"26) The FDM formulas are derived in Appendix B.III. Fig. 2.12 demonstrates an excellent agreement between the numerically simulated response curve and the analytical results. The minor graphical misfit is due to the inadequate data points in the original papers [Cooper et al., 1967; Papadopulos, 1973]. The 8 (nt = 55, tf = 1.15) value of these response curves equals to 10. This is appropriate in the inversion of the response curve of the 2D heterogeneous model, while a value of 50(nt = 100, t/ = 1.1) should be used in the inversion of of the 2D response curve from perturbed homogeneous model. Chapter 2. SIMULATION OF SLUG TEST 37 -3 - 2 - 1 0 1 Dimensionless time Figure 2.12: Simulated wellbore response (dashed lines) versus analytical solution(solid lines). From lower left to upper right, Ss decreases from 10 - 1 to IO - 1 0 . For all simulated curves, 8 = ^ = 10. Chapter 3 S Y S T E M A N D F I L T E R I N G Centering on system and filtering, we will introduce three topics in this chapter: a) The concepts of system, filtering, and the definition of a slug test conductivity filtering system; b) The properties of the slug test system filtering; c) The method to obtain the system outputs. 3.1 T h e concepts of system and filtering A system can be viewed as any process that results in the transformation of signals. Thus a system has an input signal and an output signal which is related to the input signal through the system transformation [Oppenheim et al., 1983, p35]. If the trans-formation changes the relative amplitudes of the frequency components in a signal or perhaps eliminates some frequency components entirely, the transformation is referred to as filtering [Oppenheim et al., 1983, p397]. If the filtering system is space invariant and linear, the filtering process may be represented in a convolution form, o o Y(x) = J H(x- x') X(x') dx' (3.27) —oo where H(x) is called system function or filter function, and its Fourier transform is called transfer function. The convolution is in essence a summation or superposition. A filtering system is shown is Fig. 3.1. 38 Chapter 3. SYSTEM AND FILTERING 39 INPUT SYSTEM OUTPUT X(x) H(x) Y(x) The system function epitomizes the system behavior, hence characterizes the rela-tionship between the input and the output. In a linear time/space invariant system, the system function is the system response to a Dirac input, i.e., the impulse response of the system. But in a nonlinear system, the system function can not be easily obtained, and in which case we can probe the system properties by an input of various signals and the observation of the output signals. 3.2 T h e slug test conductivity filtering system In this section we first physically define the slug test filtering system, then we discuss its appropriate mathematical representation. 3.2.1 Physical definition In Chapter 2 We have introduced the slug tests in heterogeneous core-scale(of a wellbore diameter Axw) conductivity fields, Kcore, and the numerical simulation of the wellbore response curve. If we fit the wellbore response curve by that of a homogeneous conductivity, Kalua, of the model of Cooper et al. [1967], we obtain a uniform conductivity Kslug which is thought to be somehow a representation of the heterogeneous conductivity j(core Then o n e m a y a s k ; JJ o w good is the representation? What volume of Kcore is Chapter 3. SYSTEM AND FILTERING 40 being sampled? and what is the quality of sampling(the extent of representativity) on each part of the sampled KC0Te1 More generally, what is the scale relationship between Kcore and Kalugl We utilize a spatial filtering approach to explore the scale relationship between Kcore and K"lug, by conceptualizing the scale relationship in a system and filtering framework. We regard the slug test as a conductivity filtering system which receives KCOTe input and produces Kslug output. The slug test filtering system is graphically shown in Fig. 3.13. Input Slug test filter Output Figure 3.13: Slug test filtering system receives core-scale conductivity field and produces measured-conductivity field on a larger scale. The slug test filter is what we are after. If we know the input Kcore and the output Kslug, then we can readily identify slug test filter with the well established theory of system identification, which shall be discussed in the next chapter. The input KC0Te is generated by the geostatistical simulation methods of GSLIB [37] or a spectral method elucidated in Chapter 4. The output Kalug can be obtained by two steps: First, performing slug tests at different locations in the Kcore field and recording the wellbore response curves. This is done by the numerical simulation model described in Chapter 2. Second, inverting the wellbore response curves to the conductivity values Chapter 3. SYSTEM AND FILTERING 41 Kalua. The two steps are shown in F i g . 3.14. The inversion procedure w i l l be discussed in Section 3.5. perform slug tests at different locations invert to K8"1* j£Sl l lg wellbore response Figure 3.14: Slug test scenario. The above discussion has only been concerned about the filtering of the conductivity field. It is known that specific storage is also heterogeneous in reality, thus a similar filtering may exist: S™re —• Salug. However, it is reported that specific storage is less variable than conductivity [Freeze, 1975; Dagan, 1989]. O n the core-scale, it is not unreasonable to assume S„ to be uniform. This assumption wi l l greatly facilitate the analysis of KCOTe —• Kalug filtering. W i t h both S?re and Kcore heterogeneous, the inverse method must accommodate the following difficulties: 1. The inverse estimation of Kalug and Salug can be nonunique due to the interactivity of Kalug and Sa,ug in the inversion. The heterogeneity effects of Kalug can be absorbed in Sslug and may appear pseudo-homogeneous; vice versa for Satug. It also means that the bias of one of them can propagate to the other, thus Kslug and Salug may not be representative of Kcore and S?re. Chapter 3. SYSTEM AND FILTERING 42 2. With heterogeneous S^OTe, the influence zone may be of irregular shape and differing area, in which the filtering power distribution 1 may be very erratic.2 This contra-dicts with the preferred technical requirement of the spatial filtering approach.3 For simplicity, better accuracy, and the facilitation of the technical requirement of con-ductivity filtering, we proceed with a uniform core-scale specific storage, and we examine the role of specific storage in the measurement process by varying Sg0re uniformly instead of heterogeneously. The conductivity filtering is transient process, resulting from the transient head re-covery at the wellbore and the transient propagation of the hydraulic head perturbation, Ahp(r,t), which is where the notations are the same as those in Chapter 2. At initial time t = to = 0, Ahp(r,t) = 0 for r > rw, no filtering occurs; as it starts off, the head perturbation propagates outwards from the wellbore with deceasing magnitude, H(t) \ h(r, t) \ 0, thus, the filtering area enlarges with time while the filtering power proximal to the 1 Filtering power is equivalent to filter amplitude. 2In Chapter 2 we mentioned that specific storage determines influence zone. For homogeneous S,, the influence zone is symmetric, enlarging with decreasing S,. For heterogeneous S,, one can infer the influence zone to be asymmetric and of differing area for differing S,, in which case, being convolved by the heterogeneous effects of Keore, the filtering power distribution may be drastically different for slug tests performed from one location to the other. 3The conductivity filter estimated with the spectral analysis method of Chapter 4 is an equivalent filter that can be interpreted as the one averaged from the individual filters for each slug tests sampling different core-scale conductivities; this requires the conductivity filterings of each slug tests should not be drastically different, otherwise the spectrally estimated filter may not be representative or meaningful. The different, irregular filter area and power distribution due to heterogeneous 5 £ o r e may result in fundamentally different filterings for each slug test. Thus it may reduce the representativeness of the filter identified by the spectral method. A*'(r,*)=< t > to < H(t), r = r, h(r,t), r > r, Chapter 3. SYSTEM AND FILTERING 43 wellbore decreases wi th t ime. This phenomenon shall be elaborated in the perturbation method of Chapter 5. It is pictorially shown in F i g . 3.15, where the filter functions are illustrated by three hypothetical Gaussian probability density functions 4 , / ( . , » ) - ^ e - 6 [ ' * ' ' + l*>'] (3.28) 7T Xx Xy where, Xx = Xy = 8, 16, 32, which are the characteristic filter widths. at time / , at time t. at time t. 0JJ25 Y 0.015 k Figure 3.15: Spatial conductivity filtering is a transient process. The top three show the core-scale field under filtering at different t ime t\, ti and £ 3 , where tz > ti > t\. The bottom are the corresponding filter functions. The transient filter can be obtained by inverting the response curve in different t ime intervals. The inverted Kslug in each interval lumps the spatial filtering of Kcore in that 4Here we use such a function is just for illustration convenience because of its similarity to the slug test filtering behavior, see Fig. 5.37, and explanations on page 115. Chapter 3. SYSTEM AND FILTERING 44 interval. However, this is rarely the case in engineering practice in which one observes a response in a very short time period, and one may frequently find that some parts of data record contain noise, e.g., the very early-time response. Thus, it is a usual practice to invert the response curve as a whole. This is because that the practical concern is on the gross sampling volume of a slug test in the whole duration of measurement. Considering its practical applicability, in the thesis research we adopt the method of the inversion of the whole response curve. The inversion method is equivalent to integrating the slug test filtering over the whole duration of the measurement. The identified filter function is a spatial function that represents the total sampling volume and bears different filtering power in different parts of the volume. Note that the transient nature colligates with the spatial-varying property, meaning that the spatial filtering proximal to the wellbore represents the early-time filtering, while the distant spatial filtering displays the later-time filtering. Consequently, the spatial filter illustrates a temporally integrated filtering. We will revisit the subject of temporal filtering by the perturbation approach in Chapter 5. where Yc(x), Ym(x) are the logarithm conductivities of the core-scale field Kcore and measured field Kslug, G(x) is the filter function and N(x) is the noise. The noise term accounts for numeric model misfit and random truncation error resulting from the slug test simulation method(Section 2.2, 2.3, 2.4) and the inversion routine of the wellbore response curve. 3.2.2 Mathematical representation The conductivity filtering can be mathematically represented as, (3.29) Chapter 3. SYSTEM AND FILTERING 45 The reasons for choosing logarithmic conductivity in eq. (3.29) include 1) logarithm K has a normal distribution [41]; 2) logarithm K is less variable than K; 3) the logarithm conductivity occurs naturally in the flow equation, &h dlnKdh _ 1 dh dx] + ~dx~dx~i ~ K S°W l ~ h 2 - ( 3 - 3 0 ) It can be seen that the resistivity shows up in the r.h.s of above equation. For a uniform flow in a vertically homogeneous aquifer, the effective K is a harmonic average where ^ is the basic element. The conductivity of the form, -j^ , T ^ , 5 are utilized in the literature [Gelhar, 1994; Oliver, 1990, 1992]. For the slug test filtering, it is not essential which form of K is used in eq. (3.29) as long as the form is physically meaningful and mathematically convenient. So is the expression of the filter function. One may write G{x) as G(x) = F(x)±- (3.31) x* to explore the ^ behavior reported in the pumping test literature [ Desbarats, 1994; Cardwell and Parsons, 1945]. We use eq. (3.29) in the thesis. 3.3 System inputs: core-scale fields We can probe the system filtering property by the input Kcore fields of some characteristic spatial variability (e.g., the variance, correlation, and variogram model). Because the system filter characterizes the scale relationship between the input and output, the system filter should demonstrate some features pertaining to the characteristics of input and the output. We postpone the discussion of input design to Chapter 6 until we have learned more aspects about the filtering process. 3KD is the conductivity normalized by the geometric mean. Chapter 3. SYSTEM AND FILTERING 46 3.4 Properties of System filtering In this section, we explore the system and filtering properties. 3.4.1 Stability and nonuniqueness Initial results show that the filter function in eq. (3.29) satisfies J J G(x,y) dx dy « 1, (3.32) thus, max{\gKm(x,y)} < maxjlg R~c(x, y)} < +oo (3.33) This means that the slug test system produces a bounded output from a bounded input, i.e., it is a stable system.6 The system inversion, that is inverting the system output Kslua to its input KC0Te, is nonunique. Because there are uncountable realizations of Kcore fields that can be filtered by a function G(x,y) to give an identical Kalug. 3.4.2 Space invariant property The space invariant property is similar to the time invariant property in the context of geophysical or electric engineering and it is identical to the one in image analysis, meaning that the filter function is invariant with respect to spatial translation. For the conductivity filtering, the space invariant property requires the filter volume7 and filter amplitude G(x) versus x remain the same regardless of the slug test locations in the K°°re field. Apparently, this is not strictly true. For example, Kcore filtering is limited 6Note: lg denotes base 10 logarithm in this thesis. 7Filter volume is the integral of filter amplitude in the influence zone. Chapter 3. SYSTEM AND FILTERING 47 proximal to the wellbore if the wellbore is surrounded by an annular low conductivity clay ring; while the filtering may propagate to a large volume in an isotropic sand aquifer of small variance.8 Thus, the conductivity filter can only be approximately space invariant under certain hydrogeologic conditions, e.g., isotropic media of small variance. In an absolute sense, the filtering pattern of Kcore —• f{alu9 should vary for the slug tests performed at different locations in a heterogeneous Kcore media. But, the mean or the average of the individual filters can be considered to be space invariant if the filtering function process is stationary,9 G(x-x';u>) = {G(x-x')) + g{x';u) (3.34) where G(x — x'; tn) is a realization of filter function representing the conductivity filter of slug test performed at location x'. Bracket denotes mean, (G(x — x')) is the constant mean filter and g{x'; ui) is the variation of G(x — x'; u>) at x'. The mean of G(x — x'; to) will approach the constant (G(x — x')) if we increase the number of slug test measurements, in other words, the mean of g(x'; UJ) turns to zero. If the filtering process is not stationary, (G(x — x')) may not approach a constant as more measurements are averaged. For example, in a sand-shale field [Desbarats, 1987; Bachu and Guthiell, 1990], the filtering will exhibit two distinctive patterns for slug tests performed on the sand and on the shale. The mean filter (G(x — x')) may not be meaningful. The system identification method described in Section 4.1 of Chapter 4 is developed under the space invariant assumption. The identified filter would always be an averaged filter10 from the ones of invoked by each individual slug tests. The 8In these two cases, the influence zone may be the same because it is independent of conductivity if /? = Kt/r^ is sufficiently large. What differs is the filter amplitude within the influence zone. 9The filter function can be considered to be a random process. At each location x', there is a filter function realization(labeled u), G(x — x';u). 1 0It is called an effective filter if the filtering process is stationary and the number of measurements is sufficient; or less strictly, called an equivalent filter if the filtering process is nonstationary or the number Chapter 3. SYSTEM AND FILTERING 48 representativeness of the filter depends on the extent of nonstationarity. 3 . 4 . 3 Linearity The linearity of the filtering process can be evaluated by two kinds of tests. The first is to test the scaling property or homogeneity, that is aYcX(x,y) —• aYmi(x,y), (3.35) where a is a scaling factor. The second is to test the additivity that is Ycl(x,y) + Yc2(x,y) —> Yml(x,y) + Ym2(x,y), (3.36) where Ycx(x, y), Yc2(x, y) are two different core-scale fields, and Ymi(x, y), Ym2(x, y) are the corresponding measured fields. If a process exhibits both homogeneity and additivity, then it is a linear process [Oppenheim et al., 1983, p43]. The linearity analysis of the signals K, Y = \gK, R = is shown in the Table 3.1. The subscripts c and m notate the core-scale conductivity and the measured conductivity, respectively. We use K®, Y®, R® and K®, Y®, R® denote the scaling property and additivity of variables K, Y, R. In the slug test filtering, ® is more straight forward than ®, particularly, the K® property. So we first discuss K®, then use it to analysis the other properties. Note that scaling a random variable x by a factor a changes the mean from ux to apx, variance from o~x to a2<rx; a n ( I adding two random variables x and y ends up with a mean px+y = ux + uy, and a variance o~l+y = o~x + <j* -r2crxaypxy, where pxy is the correlation coefficient. The/if® can be justified from the finite difference eq. (2.16), (2.17), (B.I.155), (B.III.171), i.e., scale Kc by a factor a is equivalent to scale Km by the same factor. It corresponds to vary the mean /xy to /iy + a. The Y® is only satisfied under the ideal condition of Yc of measurements is insufficient. Chapter 3. SYSTEM AND FILTERING 49 Table 3.1: Linearity analysis of K, IgK and j^. scaling property additivity K Y = \gK R = jt a (Kc)m = (a Kc)m lg{[(*e)m]B} ± lg{[TO°U a _ 1 {Kc\)m + {Kcl)m = (Kc\ + Kc\)m lg[(KciU + lg[(Kc2)m] I \g[(Kcl Kc2)m] I 1 I (KC)M ~ ( ^ ) M ( K c l ) m + (Kc2)m ( k e ? + J being homogeneous, by using the property of K®. The extent of Y® for heterogeneous Yc can be explored by slug tests in Yc fields with variances varying from very small (close to homogeneous) to large values. R® is satisfied thanks to the K® property. The A'®, V®, R® properties can be first appreciated in the homogeneous case in which all are satisfied. For heterogeneous fields, the extent of additivity can be explored by varying o~\ or o~Y from very small to large values. Note that Y® is satisfied if one of the fields is homogeneous, implying that if the variance of one field is small enough, regardless of the variances of the other field, Y® can be approximated. Now, we explore the extent of Y® and Y® in heterogeneous core-scale fields. A spectral method, which shall be discussed Chapter 4, page 95, is used to generate the core-scale fields, which have the spectrum S(fxify) — 2 a2 TT Xx Xy (3.37) (1 + 47T* fl Xx + ft A,)* where a2 is variance and = Ay is correlation length. For the tests of Y®, 7 fields are generated with a2 = 0.01, 0.1, 0.2, 1, 2, 4, 6, same Xx = Xy = 4Ax w and field size = 256 x 256 blocks(Axtu x Ax"). The last six fields are obtained by scaling the first field by an appropriate factor to change its variance. Slug tests are performed in every other blocks in a central area of 128 x 128 blocks in each field. Ss is chosen to be 0.01. Test results are shown in Fig. 3.16. It is evident that the extent of Y® decreases with variance. Chapter 3. SYSTEM AND FILTERING 50 For the tests of V®, 2 sets of fields are generated. Each set consists of 3 fields that have <r2 = 0.01, 0.2, 1 and the same Xx = Xy = 4Axw. The amplitude spectra between each set are identical, but the phase spectra between each set are chosen to be different. Then, respectively add up those fields of same variances between the two sets. Other parameters are the same as those of V® test fields. Test results are shown in Fig. 3.17. Clearly, we see the extent of Y® decreases with increasing variance. In other circumstances with different the specific storages and different spatial vari-ability of Yc fields, one would expect different linearity observations. 3.5 System outputs: inverse estimation of K and Sa The inverse estimation of parameters in heterogeneous media by fitting the type curves of some homogeneous models, e.g., models of Theis [1935]; Cooper and Jacob [1946]; Cooper et al. [1967], is the usual practice in various engineering disciplines. The aim is to characterize the heterogeneous properties by a few simple parameters. Even though being replete in practice, the methodology has not been adequately evaluated for the effectiveness of the inverted parameters with respect to spatial variability. Here, we only discuss the inversion of the measured response curve of the slug test in heterogeneous media. The inversion is to fit the measured response curve with that of the homogeneous model of Cooper et al. [1967]. The parameters pertaining to the fitting curve are the equivalent parameters (or so-called measured parameters) that 'somehow' epitomize heterogeneous media. The 'somehow' shall be explored in great depth in the following Chapters by the spatial filtering approach. Here we only interpret the 'somehow' from the inversion-process perspective. This is not only of great practical significance, but also is very important for the validity of the spatial filtering approach which works best for the outputs containing minimal bias. In the following, we will investigate the role of the inverse estimation variables(i.e., Chapter 3. SYSTEM AND FILTERING 51 Ym(J 4/0.01 * Yc) Yn.g 6/0.01 * Y«) Figure 3.16: Scaling property Y®: (Yc -> Ym) = • (a Yc a y m ) , where o\c = 0.01. Chapter 3. SYSTEM AND FILTERING 52 Chapter 3. SYSTEM AND FILTERING 53 conductivity and specific storage) in the process of response curve inversion, and interpret the inverse estimates with respect to the spatial variability. We introduce these issues in three step-by-step subsections: 1) The inversion objective function and its relationship to the estimation variables. 2) The estimates of two different inversion methods and their comparison. 3) Interpret the estimates of the two inversion methods with respect to the spatial variability and the role of the estimation variables in the inversion. In all of the subsections, we discuss two inversion methods which may be both appropriate: a) Estimate conductivity only, with specific storage equal to true value; b) Concurrently estimate conductivity and specific storage. 3.5.1 Objective function versus estimation variables We have discussed that the system output is obtained by fitting the measured response curve to that of the homogeneous model of Cooper et al. [1967]. The associated curve-fitting objective function, e(\gK, lgSs),n can be defined as nt e(\gK,\gSa) = 52[Hm(ti) - Ha(ti-,K,S.)]2, (3.38) i = i where Hm(ti) is the measured wellbore response curve from the simulation model dis-cussed in Chapter 2; Ha(ti; K, Ss) is the analytical solution of model of Cooper et. al. [1967], which is reproduced by the ID FDM model introduced in Chapter 2, and ti is the time with nt time steps. We find such lg K and lg SB (equal to lg K9lu9 and lg SfU9 respectively) that minimizes the objective function. This is done with a steepest decent algorithm in the MINOS [71] optimization package. From the flow equation in heterogeneous media, eq.( 2.16), we know that the wellbore response curve is a function of Kcore and Saore. Thus the response curves of various shapes 1 1 Here IgK and lgS, are utilized as estimation variables because they are less variable than K and S,. But we still use K and S, as the input parameters to the ID model of the fitting curve in eq. 3.38. Chapter 3. SYSTEM AND FILTERING 54 can be encountered in the field measurements. Then it is important to recognize that when we fit the measured curve, represented by the objective function of eq. (3.38), the estimation variables play different roles in the control of the fitting curve. Hereafter we denote the measured response curve RC and the fitting curve FC. In the inversion procedure, the changes of the estimation variables lg K and lg Ss, i.e., lg K + tfigjc and lg Ss + S\^s,? are equivalent to rescaling of K and Sa in the ID model by IO5'** and 10*'«s« respectively. From the r.h.s. of finite difference eq.(B.III.171) of the ID model, s, h? - hT'1 K At we see that scaling K will shift FC horizontally, because the scaling factor 10 5'«K can be scaled to the time, A i , i.e., \gt + ^i g AT, thus, horizontally shifting the FC of \gt versus H/HQ by S\gK- We also see that scaling S3 by 105 l« s» will change the FC shape, because K h m — hm~1 10 !«s» can be scaled to - * — ^ — , the time derivative that determines the curve shape. Thus, in conclusion, lg K determines the FC time coordinates while variable lg Sa controls the FC shape. In order to understand the relationship between the objective function and the estima-tion variables lg K and lg Sa, a hypothetical test is conducted. First, inputting a pair of K and 5 „ denoted as Ktrue, S' r u e, to the ID FDM model, yields Ha(U; Ktrue, 5^rue); second, perturbing K t r u e , Slsrue by an amount of 6 l i K , 8 h S „ gives Ha(U;Ktrue + 10*'**, u e + 10 5 l« s»); lastly calculating the objective function as e(\gK,\gS) (3.39) = £(lgi^ e + ^ , l g S r e + ^s,) nt = £ [ Ha(U;Ktrue + IO5**,Slrue + 10 s'« s») - Ha(ti-Ktrue,Slrue) ]2, SigK € [-0.15,0.15], 6lss, € [-0.15,0.15]. Chapter 3. SYSTEM AND FILTERING 55 1 A 3 I M 1.88 1.91 1.94 1.97 2 . 0 0 2 .02 2 X 6 2 . 0 8 2 . 1 1 2 . 1 4 2 . 1 0 l o g 1 0 K ^ = 2 log 1 0S s»™ = - l 1.83 1.86 1.88 1.91 I M 1.97 2 X 0 2 .02 2 -06 2 . 0 8 2 . 1 1 2 . 1 4 2 . 1 0 l og 1 0 K 1.83 1.8B 1.88 1.91 I J 4 1.97 2 . 0 0 2.02 2 .05 2 . 0 8 2 . 1 1 2 . 1 4 2 . 1 8 l o g 1 0 K ^ = 2 log10Ss»™e=-10 1 0 CO 1 0 „ I I I I ll h ll I I I I L±J -1a17 1.83 1.86 1.88 1.91 1.94 1.97 2X0 2.02 2.05 2.08 2.11 2.14 2.16 l og 1 0 K Figure 3.18: Objective function versus logarithm conductivity and specific storage: e( IgK, lgSs) versus (lgif, \gSa). Top graph, e{\gKtrue, lgS' r u e) = e(2, -1) = 0; bottom graph, e(\gKtrue, lgS' r u e) = e(2, -10) = 0. Contour level: 0.4, for e(\gK, \gSs) > 0.4; 0.0025, for e(lgK, \gSs) < 0.2. Chapter 3. SYSTEM AND FILTERING 56 Figure 3.18 shows the objective function versus logarithm conductivity and specific storage, i.e., e(\gK,\gSa) versus (IgK, \gSa). Two extreme cases of the e(lgK,lgS„) versus (lg K, lg Sa) relationship are demonstrated in Fig. 3.18, where lg Ssrue = —1, —10. In graph A, \gSaTue =—1, we see that the gradient QJZJK > d\£s. > ^ m s * n e optimum is more easily found and has a higher precision(i.e. closer to its true value) than lg Sa. In graph B, lgS^ r u e =—10, a ^ K > af£Si « 0, which means a wide range of lgS"s produces the same objective function. Any optimization method terminates the optimal searching if the objective function approaches the specified cut-off criteria. Thus, the lg Sa optima could happen to be randomly picked up from the objective function contour that has the same value as the cut-off criteria, and the lg Sa optima contain different amount of biases depending on the searching pattern of the specific optimization method used. We can also see that of graph B is greater than that of graph A, hence, with given cut-off criteria, the lg K estimation of graph B contains less bias. For intermediate specific storages, —9 < lg Sa < —2, the objective function lines will rotate approximately from NW 30°(graph A) to N90°(graph B). In summary, with a given cut-off criteria, we observed: 1) lg K optimum is more easily found and has a higher precision than lg Sa optimum. 2) For smaller specific storages, the lg K contains less bias whereas lg Sa estimates is more strongly biased. Note that these observations are based on the fitting to the response curve of a homogeneous aquifer, and they are not necessarily true for the fittings to those response curves of heterogeneous aquifers. For example, for the fit to a response curve of some particular shape, df*s may be greater than d ^ K so that lg Sa is more easily found and has a better precision. The fit to the response curves of heterogeneous aquifers can be easily carried out by a simple numeric routine, but the interpretation of Kalu9 and SfU9 in terms of the spatial heterogeneity is not straightforward, and it is very complex to quantify the role of the estimation variables in the inversion process. In the following, in Section 3.5.2 we first Chapter 3. SYSTEM AND FILTERING 57 present the estimates of two different inversion methods and their comparisons, then in Section 3.5.3 we interpret the estimates of the two inversion methods with respect to the spatial variability and the role of the estimation variables in the inversion. 3.5.2 Comparison of the estimates of different inversion methods Here, we introduce two inversion methods denoted as a and b: a) Estimate lg K only with lg Sa = lg Slrue(hereafter we denote the estimate lg if as lg Ka). b) Concurrently estimate IgK and lgSs(hereafter we denote the estimates as \gKh and lgS*). We know that the RC measured in a heterogeneous conductivity field usually does not exactly conform to the FC of an equivalent homogeneous media. Thus, both the estimates of the two methods are more or less biased. For method a, the e(lg K) of the shape-misfit between Hm(ti) and Ha(ti; K) can not be reduced, whereas for method b the shape-misfit can be minimized by adjusting the shape-control factor lg Sb. However, the lg S^ deviating from lg Slrue is a biased estimate, and the bias may pollute the \gKb with larger bias than lg Ka due to the nonuniqueness (or interactivity) of lg Kb and lg Sb in the inversion. We perform slug tests and measure the RCs in an isotropic heterogeneous Kcore field with an exponential variogram of crfsK = 0.25 and correlation length of A = 4 Ax™(wellbore diameter), 7 W = "J* • [ 1 - ^ *>]• (3-40) The field is generated by a sequential Gaussian simulation method of GSLIB [37], and it comprises 256 x 256 square blocks. A number of Sa values (true Sa) are assigned to the field, Ss = 10"1, 0.05, I O - 2 , 0.005, 10"3, 10~4, 10"5. Slug tests are performed in every blocks in an central area of 64 x 64 blocks. Figure 3.19 shows the comparison of the histograms of the lg Ka and the ones of lg Kb and lgS*, with true S3 varying from 10 - 5 to 10_ 1. We see that: 1) Both \gKa and \gKb shrink towards the mean of the core-scale field, pl KCOre « 0. 2) The shrinking is skewed, Chapter 3. SYSTEM AND FILTERING 58 with the high K end shrinking more strongly than the low K end, and lg Kb being more skewed than lg Ka. The skewness becomes more pronounced with decreasing true Ss. 3) The distribution of lg S% is not normal to the true S3, spreading out from the true Ss and asymmetrically skewing towards the low value. Larger spreading occurs when the true Ss decreases. Fig. 3.20 displays the cross plots lg Sb versus lg Kb, with true Sa varying from 10~5 to 10 - 1. We see that: 1) Most lgS* estimates are biased,12 and the magnitude of bias increases with decreasing true Sa. 2) For decreasing true Ss, the high conductivities shrink more strongly towards the 0, but, the low conductivities are somehow persistent. Figure 3.21 illustrates the comparisons of the lg Ka with the lg Kb by the cross plots lgiY" versus \gKb. We can roughly see that: 1) For true Ss = 10 _ 1, \gKa is mostly smaller than lgKb, but for Ss < 10_ 1 the high value of \gKa is greater than that of lgKb. 2) For decreasing true Ss, the high conductivities shrink more strongly towards the mean p l g K C O r e ^ 0, but, the low conductivities are somehow persistent; and the high conductivities in lg Kb is more strongly reduced to the p l g K C O r e ~ 0 than is the lg Ka. Fig 3.22 shows the cross-plots of [(\gKa) - {\gKb)\ versus [(\gSbs) - (lgS'rue)]. We see that: 1) When \gSba equals to lgSf u e , \gKb collapses to \gKa. 2) Whenever \gSbs < lgS' r t i e, \gKb > \gKa. 3) The derivative of [(lgSbs - lg^ r u e)] with respect to [(lgKa) -(lg/f6)] increases with decreasing l g ^ " 6 , which means a given amount of bias of IgS^ will cause less bias in the lg Kb estimation with decreasing lg 5^rue. In summary, these four types of graphs showed two major points: 1) For decreasing Ss, lg Ka and lg Kb skewly shrink towards // l g K C O r e w 0, while lg Sb asymmetrically spreading out from the true Ss. The low-value conductivity is more persistent than the high-value conductivity. 2) Whenever \gSbs < lgSj r u e , lg Kb >\gKa. 1 2The bias refers to |lg5* - lg5jrue|. Chapter 3. SYSTEM AND FILTERING 59 Figure 3.19: Comparison of histograms of the lg K estimated with true lg Sa and the ones of the IgK concurrently estimated with \gSs. Chapter 3. SYSTEM AND FILTERING 60 Ym* versus loglO(Ss) hYm*J: from concurrent K and Ss inversion -« -10 : 1 1 1 : , .V*3Blg?F'-. • Ss = 10-2^ i . I Ss = 10"1 Figure 3.20: Cross-plots of lg K concurrently estimated with lg Ss versus estimated lg S9 The dashed lines represent the true lg Ss. Chapter 3. SYSTEM AND FILTERING 61 Y m a versus Ym* |Ym ap from inverting K alone from concurrent K and Ss inversion Figure 3.21: Comparison of IgK estimated with true lgSs with \gK concurrently esti-mated with lg Ss. The inclined lines are of unit slope. Chapter 3. SYSTEM AND FILTERING 62 (IgK estimated with true lgSs) - (IgK estimated concurrenUy with lgSs) Figure 3.22: Cross-plots of [(IgA' estimated with true lg5«) - i}gK estimated concur-rently with lg5s)] versus [(estimated lgSs) - (true lg5s)]. 3.5.3 Interpretation of the estimates of different inversion methods In order to explain the causes of the skewed histograms and the differences between the two types of inverse estimates, first in Section 3.5.3.A we discuss some extreme values of the inverse estimates of the measured data. Then in Section 3.5.3.B we investigate the effects of spatial conductivity variability and the effects of specific storage on the inverse estimates. Lastly, in Section 3.5.3.C we explore the influence of low and high conductivity at different locations from the wellbore. A summary is given at the end of Chapter 3. SYSTEM AND FILTERING 63 Section 3.5.3.C. 3 . 5 . 3 . A Extreme values of conductivity estimated with true Ss We discuss the maximum and minimum of the conductivity estimates inverted with method a. We denote the core-scale conductivities corresponding to these extrema as Khi9h and Klow respectively. The RCs(response curves) of these Khigh and Klow are also inverted with method b. The skewness, K™, and the differences between the estimates of the two inverse methods, AK and ASa, are defined in table 3.2. The skewness is a parameter used to describe the extent of skewing of the histogram of the measured fields with respect to their means. Table 3.2 shows the estimates from the RCs simulated Table 3.2: The maximum and minimum estimates, skewness K*, AK and ASS. K and Sa of fitting curve * varying K only varying K and Sa Cjtrue K K" * K Ss K™ AK t ASa * J^high -1.0 -2.0 -3.0 -4.0 -5.0 1.0208 -0.5757 0.6773 -0.6855 0.3724 -0.7954 0.2042 -0.8264 0.1467 -0.7867 0.6863 -0.3151 -0.0553 0.3115 -0.8557 -0.5349 -0.1411 -0.9315 -1.1200 0.0035 -2.5212 -0.6634 0.0074 -3.6017 -0.6766 0.3345 0.6849 0.3658 1.1443 0.5135 2.0685 0.2007 1.4788 0.1393 1.3983 -1.0 -2.0 -3.0 -4.0 -5.0 -1.7143 -1.4806 -1.2856 -1.1484 -1.0512 -0.8594 -3.8528 -0.9642 -4.7528 -0.9789 -5.7589 -0.7847 -9.0000 -0.8018 -9.2789 -0.8549 -2.8528 -0.5164 -2.7528 -0.3067 -2.7589 -0.3637 -5.0000 -0.2494 -4.2789 *: all are logarithmic data; •: K" = Khi9h + K'ow - 2 uKCOre, for Khi9h > 0 > fiKCore = -0.0589; = Khi»h + Klow, for Khi9h < 0 < pKCOre; pKCOre = mean of core-scale field; f: AK = (K of varying K only) - (K of varying K and S„); t- ASS = (S. of varying K and S.) - (S]™). Chapter 3. SYSTEM AND FILTERING 64 with total time tmax = P = 10, where Kn equals the smaller of the geometrical mean of the core-scale conductivities at the inner zone of the simulation model, and the geometrical mean of the conductivities at the four blocks around wellbore. Such a tmax is able to catch the whole recovery in Kh%9h and Klow as well. A negative K* means the histogram skews to values less than the mean. We see that all the estimates are negatively skewed, the extent of skewness is proportional to decreasing true Sa. AK and ASa represent the differences of the inversion methods a and 6. They are seen to be positively correlated. Graphic illustration of K* versus Sa and AK versus ASa shall be given in the next subsection where a series of comparisons are displayed. We graphically show, in the left column of Fig. 3.23, the RCs and the FCs(fitting curves) of Kh%ah and Klow by varying K only; in the right column, the RCs against the FCs by varying both K and S„. We see: 1) For Klow, FC of method a always crosses the RC at approximately H/H0 = 0.5, and FC is lower than RC at H/H0 < 0.5 while high than RC at H/H0 > 0.5; The FC of method b fits the RC quite well. 2) For Khi«h, the relative positioning pattern of the FC and RC in method a is opposite to that of Klow except the very early-time FC of Ss = 0.1. FC of method b does not closely fit RC for S3 > 0.001. Chapter 3. SYSTEM AND FILTERING 65 continued on next page • Chapter 3. SYSTEM AND FILTERING 66 solid line: simulated response curve solid line: simulated response curve dashed line: fit with varying K only dashed line: fit with varying K and Ss LoglO(t) LoglO(t) continued on next page • Chapter 3. SYSTEM AND FILTERING 67 solid line: simulated response curve solid line: simulated response curve dashed line: fit with varying K only dashed line: fit with varying K and Ss LogltXt) LoglO(t) Figure 3.23: The solid lines are the simulated response curves. The dashed lines are the fitting response curves. Left column represents the fit by varying K only. Right column is the fit by varying both K and Ss. Chapter 3. SYSTEM AND FILTERING 68 3.5.3.B Effects of spatial conductivity variability with different specific stor-ages on the inverse estimates In this section, we investigate the K™ versus Sa and AK versus ASa relationships ob-served in the inversion of Kh%9h and Klow through a series of tests in conductivity fields of hypothetical variability. Because the KCOTe (core-scale conductivity) proximal to the wellbore has the strongest influence on Kalug (measured conductivity from slug test), the near-wellbore Kcore corresponding to the maximum Kalu9 may be the highest Kcore conductivities; and due to the; existence of correlation in Kcore, one may expect the Kcore somehow varies, away from the wellbore, from the near-wellbore high conductivities to some low conductivities, and vice versa for the Kcore associated with the minimum Kalu9. Consequently, we design a hypothetical radially symmetric spatial variability varying, according to some function, from lg Kcore(r = rw) = ± 2 to \gKcore(r = rmax(Sa)) = 0, where rmax(Sa) is the influence radius dependent on Sa. This variability should cover the practical ranges involved in the previous sections. The following is the designed continuous functions, The a' and c' are symmetric about the middle point of b'. The nice property of these functions is the exponent ct which is very flexible to adjust the variability. A series of discretized functions can be obtained if we use i of r(i) as the varying variable, the above becomes, max w (3.41) (3.42) (3.43) (3.44) Chapter 3. SYSTEM AND FILTERING 69 c: Yc(Ar(i)) = b: Yc(Ar(i)) = (3.46) (3.45) i = l , . . . , n, where Ar(i) represents the ring r(i) — r(i — 1), and n is number of rings. If we set n = 100, the function exhibits another kind of variability. Functions a, b, c, d are used in the tests, which are shown in top two graphs in Fig. 3.24. Test a represents a lg K slowly varying from the wellbore, test c a IgK fast varying from the wellbore, b a linear log variation, d represents a exponent linear variation, i.e., K varies linearly. In equations (3.44) to (3.47), we denote the positive and negative Yc (Ar(i)) as K> and KK, respectively. We perform slug tests in K> and K< with different 5", values, see the upper-left graph of Fig. 3.24. For each Sa, the summation of the measured conductivity in K> and that in K< is what we called skewness K™. A non-zero KM is a measure of non-linearity(non-additivity), in which case the measured conductivity of K> and K< are not symmetric about the mean. We conduct slug tests and the analyses of K* for all of the hypothetical variability cases a, 6, c, d and for S3 = 10 _ 1, 10~2, 10 - 3, 10~4, 10 - 5. For results comparability, we choose the r m a x of Ss = 10~5 and an empirical grid spacing of Ar(i) for all the Ss cases. The conductivity in the ring Ar(i) is given by the above formulas. The bottom two graphs of Fig. 3.24 show the skewness K™ versus Ss of these tests, compared with the K** of the measured data observed in Section 3.5.3.A, for both inver-sion method a and b. The skewness represents the extent of additivity Y® with respect to different conductivity variations and specific storage values. For a perfect Y®, we (3.47) Chapter 3. SYSTEM AND FILTERING 70 should observe K™ = 0. In method a, we see: 1) All K* < 0, meaning that Y® is not satisfied. 2)The linear-variation test 6 is the closest to the measured data. The lg K of test a varies too slowly to mimic the measured case, while c and d too fast. 3) The K* of the slowly-varying test a is the smallest for Sa > 10 - 4. This is consistent with the rule that the extent of Y® is proportional to decreasing variance. 4) The KM of the linear exponent test d does not change with Sa, it may be due to the fact that the fast-varying portion of d is influenced by the slug tests of all the Sa cases, and it may follow a long slowly-varying portion in which the extent of Y® is similar, i.e., K* is similar for all Sa cases. Fig. 3.25 shows AK versus ASa. All the graphs demonstrate the positive correlation. Test b is the closest to the measured data case, which is consistent with that in our observation in K™ versus Ss. Fig 3.26 shows the FC and RC of test d. Three cases are considered: 1) K>: K1 > K2 > • • • > Kn = K°; 2) K°: K1 = K2 = • • • = Kn = K°; 3) K<: K1 < K2 < • < Kn = K°. We see: 1) For K<, FC of method a always crosses the response curve(RC) at ap-proximately H/HQ = 0.5, and FC is lower than RC at H/H0 < 0.5 while high than RC at H/Ho > 0.5; FC of method b fits RC quite well. 2) In method a, the relative positioning of FC and RC of Ky is opposite to that of KK. 3) The skewness K* is displayed by the difference (horizontal distance between the FC of K> and K°) — (horizontal distance between the FC of K< and K°). The 1) and 2) are consistent with our observation in the FC and RC of the measured data case, which means the designed test is reasonable. The FC and RC of tests a, b, c display the same pattern. The results of tests a, b, c and d are listed in Table 3.3.. Chapter 3. SYSTEM AND FILTERING 71 loglO(K') versus i k ) g l O ( K ' ) = I curve a I + - l o g l O C K 1 ) { 1 - [ ( i - l ) / ( n - l ) ] « } ; I curve b I + - l o g l O d c 1 ) [ 1 - ( i - l ) / ( n - l ) ] ; I curve c i + - l o g l O ( K . l ) { 1 - [ ( i - l ) / ( n - l ) ] a } ; I curve d~\ + — k > g l O ( i / n ) ; l o g l O C K 1 ) = 2 ; a = 5; i = 2 , n ; n = 100 -3 I >-10 20 30 40 50 60 70 80 90 100 invert K alone: (K^= K* + K* ) versus Ss invert K and Ss: (K^= K* + K=* ) versus Ss Ss Ss Figure 3.24: Conductivity variation functions and skewness of tests a, b, c, d. Chapter 3. SYSTEM AND FILTERING 72 (logKfl - logK*) ^ (logSs* - logSs*"') -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -0.8 -0.6 -0.4 -0.2 0 0.2 Figure 3.25: AK ~ ASS of conductivity variation functions a, b, c, d. 0.4 0.6 Chapter 3. SYSTEM AND FILTERING 73 Case UK*)'. K i > K 2 > ••• > K n = K° loglO(Ki) = loglO(K«) - loglO(i/n) loglO(K°) = 0, n = 100 Case2(K°): K i = K 2 = - = K n = K° loglO(K>) =logl0(K») logl0(K°) = 0 ,n = 100 CoseJtK*): K1 <K 2<~<K" = K° loglO(K') = Iogl0(i/n) + loglO(K«) loglO(K») = 0, n = 100 solid line: simulated response curve dashed line: fit with varying K only solid line: simulated response curve dashed line: fit with varying K andSs continued on next page Chapter 3. SYSTEM AND FILTERING 74 solid line: simulated response curve dashed line: fit with varying K only solid line: simulated response curve dashed line: fit with varying K and Ss 4 -3 LogHXt) continued on next page Chapter 3. SYSTEM AND FILTERING 75 solid line: simulated response curve solid line: simulated response curve dashed line: fit with varying K only dashed line: fit with varying K and Ss LoglO(t) LoglOft) Figure 3.26: Fitting and response curves of Ky, KK and K° of test d. Chapter 3. SYSTEM AND FILTERING 76 K a n d S, o f f i t t i n g cu rve + varying K only varying K and S, gtrue K K * K S, K" AK T A S , * -1.0 1.9860 -0.0040 1.9510 -0.9003 -0.0150 0.0350 0.0997 a -2.0 1.9520 -0.0090 1.8570 -1.5860 -0.0327 0.0950 0.4140 K> -3.0 1.8800 -0.0240 1.6652 -1.8113 0.0002 0.2148 1.1887 -4.0 1.7490 -0.0773 1.3499 -1.5917 -0.3145 0.3991 2.4083 -5.0 1.5250 -0.2180 0.8710 -1.1209 -0.5850 0.6540 3.8791 -1.0 -1.9900 -1.9660 -1.0610 -0.0240 -0.0610 a -2.0 -1.9610 -1.8897 -2.3432 -0.0713 -0.3432 K< -3.0 -1.9040 -1.6650 -5.2850 -0.2390 -2.2850 -4.0 -1.8263 -1.6644 -6.2745 -0.1619 -2.2745 -5.0 -1.7430 -1.4560 -10.0000 -0.2870 -5.0000 -1.0 1.4040 -0.1010 0.9524 -0.0900 -0.2276 0.4516 0.9100 b -2.0 1.1830 -0.1470 0.8056 -0.8021 -0.1032 0.3774 1.1979 K> -3.0 0.9848 -0.2142 0.6214 -1.3235 -0.2961 0.3634 1.6765 -4.0 0.8098 -0.2882 0.4353 -1.7155 -0.4954 0.3745 2.2845 -5.0 0.6525 -0.3615 0.2607 -2.0357 -0.4665 0.3918 2.9643 -1.0 -1.5050 -1.1800 -2.1250 -0.3250 -1.1250 b -2.0 -1.3300 -0.9088 -5.3320 -0.4212 -3.3320 K< -3.0 -1.1990 -0.9175 -6.1270 -0.2815 -3.1270 -4.0 -1.0980 -0.9307 -6.6302 -0.1673 -2.6302 -5.0 -1.0140 -0.7272 -10.0000 -0.2868 -5.0000 -1.0 0.4353 -0.2941 -0.0948 -0.0067 -0.2598 0.5301 0.9933 c -2.0 0.2332 -0.2810 0.0015 -1.1830 -0.0966 0.2317 0.8170 K> -3.0 0.1443 -0.2514 0.0054 -2.1940 -0.1317 0.1389 0.8060 -4.0 0.1014 -0.2218 0.0021 -3.1750 -0.1768 0.0993 0.8250 -5.0 0.0782 -0.1952 0.0123 -4.2750 0.0258 0.0659 0.7250 -1.0 -0.7294 -0.1650 -3.6230 -0.5644 -2.6230 c -2.0 -0.5142 -0.0981 -5.2220 -0.4161 -3.2220 K< -3.0 -0.3957 -0.1371 -5.5597 -0.2586 -2.5597 -4.0 -0.3232 -0.1789 -5.7373 -0.1443 -1.7373 -5.0 -0.2734 0.0135 -9.9940 -0.2869 -4.9940 -1 0.5734 -0.1100 0.2500 -0.3312 -0.0985 0.3235 0.6688 d -2 0.4292 -0.1042 0.2311 -1.2805 -0.0326 0.1981 0.7195 K> -3 0.3371 -0.1050 0.1718 -2.0693 -0.0349 0.1653 0.9307 -4 0.2721 -0.1067 0.1126 -2.7677 -0.1265 0.1596 1.2323 -5 0.3406 -0.1067 0.1797 -3.4237 0.0093 0.1609 1.5763 -1 -0.6834 -0.3485 -2.1481 -0.3350 -1.1481 d -2 -0.5334 -0.2637 -3.6869 -0.2697 -1.6869 K> -3 -0.4421 -0.2067 -5.2189 -0.2354 -2.2189 -4 -0.3789 -0.2391 -5.6513 -0.1398 -1.6513 -5 -0.4470 -0.1704 -9.7245 -0.2766 -4.7245 *: al are logarithmic data; •: K " = K> + K<; t: AK = (K of varying K only) - (K of varying K and 5,); J: AS, = (K of varying K and S,) - (S*srue )• Table 3.3: K", AK ~ ASS of tests a, b, c, d. Chapter 3. SYSTEM AND FILTERING 77 3.5.3.C Influence of low and high conductivity at different locations from the wellbore on the inverse estimates In this section, we discuss the influence of low and high conductivity at different locations from the wellbore on the inverse estimates. In other words, we shall explore the KM versus Sa and AK versus ASa relationships in space. We employ a numeric perturbation idea, shown in Fig. 3.27, perturbing the homogeneous conductivity Yc(x) by an amount of is carried out, one block per time, on all blocks on one centerline, using the FDM in-troduced in Chapter 2. We fix the true Sa = 0.1 and observe the spatial relationships K™ versus rD, AK versus rD, and ASa versus rD, where rD = r/rw is the dimensionless perturbation location. For Sa equal to other Sa values, similar observations should be expected. Several typical RCs of /ifH a n d near the wellbore or distant from wellbore are shown in Fig. 3.27. The positioning of FC to RC, of # H a n d K^\ for both inversion method a and b, should exhibit the same pattern as we observed in the Klow and Kht9h of measured data case or the K< and K> of the test cases. Compared with the type curves of the homogeneous RC, we expect, for inversion method b, Ss shall be underestimated in the FC of near-wellbore low K, while overestimated in the FC of distant low K, and we should observe the opposite in the high K case. The location that Sa is correctly estimated is where such a perturbation has identical influence on the early- and later-time RC, i.e, x f2p, x e o p 0 < a < 1 (3.48) where fl p is the area of one grid block. We use Kt ' and K^ to notate negative perturba-tion Yc(x) — b~Yc(x) and positive perturbation YC(X) + 6YC(X), respectively. The perturbation Chapter 3. SYSTEM AND FILTERING 78 it doesn't alter the shape of the RC of the homogeneous case, only horizontally shifts it. The perturbation is independent of Yc according to the K® property (Yc(x) ± 8yc(x) is equivalent to 10±5y<:(*>.Kc), thus we can assume Yc = 0 without losing generality. We choose a =0.99, 0.95, 0.9, 0.8, 0.7, 0.6, 0.5 and S, = 0.1. The results of the tests K[~] and are shown in Fig. 3.28 and Fig. 3.29. These plots are the inverse estimates versus the perturbation location. We see: 1) Shown in the first 4 graphs, K™ decreases with a, and K™ « 0, for a > 0.9. It means that K® is approached if the changes of conductivity by the perturbation are less than 10% of the homogeneous conductivity. This somehow quantifies the previous linearity analysis that K® is approximated if the variance is sufficiently small (close to homogeneous). 2) The comparison of the results of inversion method a and b are shown in the last two graphs. The graph AK(K^) versus rD shows versus r„. In comparison of the graph ASS(K^) versus rD which is [Ss(KM)b - tvneSs] versus rD, it is clear that AK(K^) and ASa(K^) are positively correlated. 3) At the only location, rD = 1.8, ASa(K^) = 0, i.e., Sa(K^)b is correctly estimated, which is where the perturbation has identical influence on the early- and later-time RC, thus doesn't alter the shape of the RC of the homogeneous case. At this location, AK(K^) = 0, i.e., K{K^Y = K{K^)b. For other true Sa values, one may expect that the rD where Sa(K^)b is correctly estimated shall increases with decreasing true Sa. To summarize the hypothetical tests in the preceding two sections, in Section 3.5.3.B we observed: 1. versus Sa of method a is negatively skewed, and the extent of skewness of the versus Sa with the linear log conductivity variation is proportional to decreasing true Sa. versus Sa of method b exhibits erratic changes which result from the noises of lg Kb due to the interactivity lg Sb versus lg Kb in the inverse estimation. Chapter 3. SYSTEM AND FILTERING 79 2. AK versus ASa conforms a positive correlation. In Section 3.5.3.C we observed: 1. K* versus rD decreases with perturbation magnitude, and K* versus rD approaches zero if the perturbation about is less than 10% of the unperturbed conductivity. 2. AK(K^) versus rD positively correlates with ASa(K^) versus rD. The explanation for the first observations is that the slug test filtering system exhibits nonlinearity. The extent of nonlinearity is proportional to the variance and the decreasing Sa of the heterogeneous aquifer. The reason for the second observations is the interactivity of lg Sb„ versus lg Kb in the inverse estimation in which lg Sb controls the shape of FC in the fitting to the RC which can be of various shapes due to the heterogeneity. In rest of this thesis, we choose the inversion method a to estimate Ym. Chapter 3. SYSTEM AND FILTERING 80 homogeneous response curves -3 - 2 - 1 0 1 LoglOCrt/r_w*2) heterogeneous response curves heterogeneous response curves high K lowK T 1 1 1 1 I r LoglO(t) LoglO(t) Figure 3.27: Perturbation set-up, and the repsonse curves of low/high K perturtbations near wellbore and distant from wellbore. Chapter 3. SYSTEM AND FILTERING 81 Measured conductivity and specific storage in homogeneous field with a small area of perturbation versus dimensionless perturbation location r D = r/rw Notations: 1) solid lines: negative perturbation K[-] = log(K) - k>g(l/a) 2) dashed lines: positive perturbation K[+] = log(K) + log(l/a) Graph labels: from the zero line to its two sides, a = 0.99,0.95,0.9,0.8,0.7,0.5. @ invert K alone l@invert K and Ss K I _ ] - r„, and K M ,» R „ K W ~ r D , and KW ~r„ 4.0e-04 2.0e-04 \-0.0e+O0 r -2.0e-04 h -4.0e-04 -6.0e-04 4.0e-04 Z0e-04 0.0e+00 -Z0e-04 -4.0e-04 -6.0e-04 1 (tf*=KH + K r + ] )~r D H 4.0e-O4 2.0e-04 0.0e+O0 -2.0e-04 r H -4.0e-04 r--6.06-04 u 1 4.0e-04 h 2.0e-04 r 0.0e+00 \-H -2.0e-04 r H -4.0e-O4 h -6.0e-04 l - L ( l ^ = K H + K W ) ^ r D Figure 3.28: ~ rD, K™ ~ rD, where rD = r/rw, of inversion scheme a and b. True & = 0.1. Chapter 3. SYSTEM AND FILTERING Figure 3.29: A i f * ~ r c , ASS ~ rD, where rD = r/rw, of inversion scheme a and b. True 5S = 0.1. Chapter 4 S Y S T E M I D E N T I F I C A T I O N : S P E C T R A L A N A L Y S I S Spectral analysis is known as the estimation of the power spectrum of a stationary random process, and it is also applicable to deterministic as well as nonstationary process with appropriate modifications [Mohanty, 1986; Priestly, 1981]. Using spectral estimates, in Section 4.1 we identify the system filter by a Wiener filtering method. The Wiener filtering method requires the estimation of power spectrum, so we discuss the spectrum estimation methods in Section 4.2. In the last section, we test the Wiener filtering approach by the estimation of a number of hypothetical filters. 4.1 Nonparametric filter function identification: Wiener filtering We have conceptualized in Section 3.2.2 that the spatial filter can be represented as Ym(x) = J G(x- x') Yc(x') dx' + N(x) (4.49) where where Yc(x), Ym(x) are the logarithm conductivities of the core-scale field and measured field, G(x) is the filter function and N(x) is the noise. In this section, we use a Wiener filtering approach to deconvolve the system filter G(x). Taking Fourier transform of (4.49) and using hat to indicate Fourier transformation, we get Ym = GYC + N. (4.50) Let's define, Ym = GYC = Ym - N, (4.51) 83 Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 84 G = YmQ (4.52) where F m , G are optimal estimates of Ym, G, and Q is an optimal filter that minimizes 2 G - G which is G-G YmQ N, Ym are assumed uncorrelated so that the previous expression reduces to 1 Yc Yr Q - l + N \Q 2 I -(4.53) Differentiating equation (4.53) with respect to Q, and setting the result equal to zero gives Q = Ym 2 Ym 2 + —12 N\ (4.54) This is the formula for the optimal filter Q. Thus the optimally estimated transfer function is YmQ G + Yc N (4.55) Multiplying byFc*, dividing by Ym , and substituting with \YC G , yields Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 85 We rewrite above formulation as (4.56) where S'cm(/) is the cross spectral density of Yc and Ym, Scc(f), SNN{/), SGG(I) ARE spectral densities of Yc, N, G, respectively, and f is spatial frequency. Eq. 4.56 can also be deduced from the similar work of Sondhi [1972] and Helstrom [1967] in image analyses. The Wiener filtering is performed by an iterative procedure, shown in Figure 4.30. 4.2 Spectrum estimation The spectral estimate of minimum bias and lowest variance is the goal of various spectrum estimation methods. Centering on bias and variance, in this section we first discuss the sources of bias, and the dilemma existing in the spectrum estimation; secondly, we discuss the 'traditional' periodogram and data tapering method; lastly, a multitaper method, which is an improved version of the periodogram and data tapering method, is presented and verified through some hypothetical filtering tests. 4.2.1 Bias and di lemma of spectrum estimation We define the bias, Bs(f), of a spectrum estimate 5jv(/)5 as [Mohanty, 1986], Bs(f) = E{SN(f)-SN(f)}, (4.57) where «SV(/) is the true spectrum. Understanding the sources of biases in a spectrum estimation procedure by the spectral method will be of great help to minimize the bias and improve the resolution or accuracy of the spectral estimation. The typical steps in a spectrum estimation procedure and their effects on the spectral estimates are illustrated in Fig. 4.31. We see that bias can be introduced in every stage 6{f) -Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 86 Ye spectrum Cross-spectrum ScmU) = ?e(f) ?Mf) Initial signal/noise ratio fgj^J = initial value I ^ Transfer function Ql f\ — 5<=-(/) Filter spectrum SGG{I) = G(f) G{f) I Estimated Ym(f) Ym(f) = hf) Ym(f) E^(/)(,+1) 1=0 I I ^Noise N(f) = Ym(f) - Ym(f) r _E^ 0 G(/ ) ( ' + 1 >-E;* 0 G(/ ) ( ' E ^ o ^ ( / ) ( i + 1 ) I Noise spectrum SNN(I) = N(f) N*(f) NO Figure 4.30: Filter identification by an iterative Wiener filtering procedure. Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 87 Data record of finite length Data tapering D F T Averaging by spectral window XN(n) q(n) YN(n) I YN(n) • ( limited resolution, limited reliability, end effect. I YN(n) _ J _ / o * M / ) U ( / - / ' ) d/'(implicit) J A T (/) A(/ - / ' ) (taper) T Smoothing P S D - Mi sharpen spectral window, reduce leakage, increase variance. high frequency aliasing, end effect. reduce leakage, increase variance. reduce variance, decrease resolution. Figure 4.31: Spectrum estimation procedure by spectral method. in the spectrum estimation procedure, which typically are the data tapering in the lag-window, discrete Fourier transformation(DFT), spectrum averaging in spectral-window, and spectrum smoothing. Data tapering discards part of data record and unevenly weights the data in the lag window. Aliasing occurs when the Nyquist frequency, is smaller than the cut-off frequency, / c , in the original signal, i.e., 2 A < U The implicit DFT periodicity folds the spectral energy at / > ^ back to the energy at [0, 2 ^ ] . Spectrum averaging by the spectral window function, the DFT of the lag window, causes uneven averaging in the main lobe of the spectral window, while in the small side lobes it introduces leakage from spectral peaks. Spectrum smoothing is performed Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 88 through a moving-window in which adjacent estimates are arithmetically averaged. The quality of an estimated spectrum can be judged through the evaluation of esti-mation bias and spectrum variance. The smaller bias and the lower variance, the better the PSD estimate. However, these two goals usually can not be achieved simultaneously by a single data taper method. This can be seen from the fundamental result [N ewland, 1993, P108], a 1 m y/B7N' (4.58) 1 /-oo Be = ^ — , and / A(f) df = 1, (4.59) where o, m, Be, N are standard deviation, mean, effective spectral-window width and the length of data record. A(f) is the spectral window function. Eq.(4.58) holds valid if spectrum changes slowly over frequency intervals of order 1/N, i.e. provided that the record length is long enough to resolve adjacent spectral peaks. Eq.(4.58) also applies to the measurement made with an analogue spectral analyzer and applies approximately to the measurement of the autocorrelation function. From Eq.(4.58), we can now appreciate the common dilemma that is, for a high resolution spectrum(i.e. small bias), Be must be small, while for good statistical relia-bility(i.e. low variance) Be must be large compared with jj, the Rayleigh frequency. For an example, the Blackman-Tukey lag window [Mohanty, 1986], , \ ( l+cos(^) ), \k\<M, X (k) = { 2 V II— ( 4 6 0 ) 0, otherwise. It gives a bias of [Priestley, 1981], Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 89 and a variance as ' M J S(w)2, w^O,±ir, \im[N Var(S(w))] = ( 4 " 6 2 ) M f S(w) , otherwise, where w is angular frequency and S(w) is the power spectrum density. We can see that when M is large(2?e is small), the bias is small and the variance is large. If M is small(5e is large), the bias is large and the variance is small. In general, a single-taper methods are plagued by the trade-off between bias and variance, however, multitaper method provides "tricks-of-the-trade" and yield robust spectrum estimates, which we shall discuss after the next section. 4.2.2 Periodogram and data tapering The periodogram can be defined as [Mohanty, 1986], W) = jj |^v(/)| 2 (4.63) where XN(/) is the DFT of data series Xjv(n). The periodogram is an asymptotically unbiased estimate of the true power spectrum density, SN(/), in the sense that [Mohanty, 1986], lim E{IN(f)} = SN(f), (4-64) but, it is an inconsistent PSD estimate because lim Var{IN(f)} ^ 0, (4.65) which means the spectrum variance does not drop off as we increase the number of measurements. Although the periodogram does not have any data tapering bias, i.e., it uses the exact whole data record, the periodogram can be shown to be equal to the true PSD Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 90 S(f) convolved with an implicit spectral window that is the DFT of a boxcar function Iljv(n) of a unit height and a length of N, Iljv(n) = 1, n = 0,1,2, N - 1 fi = su^N) = JV_rinM x = x / J V 7 T / X S(f) = jH |n(/-/ ') | 2 S(f')df (4.67) where H(/) is the sine function and S(f) is the estimated PSD. Eq.(4.67) shows that H(/) unevenly weights S(f) and causes leakage to the roll-off frequency. Therefore we can use a tapered data window to smooth the data at each end of the record before it is analyzed. This has the effect of sharpening the spectral window, hence reducing leakage. The tapers, Hanning, cosine, Parzen, etc., belong to this kind. The PSD estimate of the tapered data relates to periodogram, S{f) = r A(f-f')\2 W)df (4-68) I00 A(f-f')df = 1, (4.69) J — CO where A(f) is the DFT of the data taper a(n). However, as long as only a single taper is used, there will be a trade-off between the resistance to spectral leakage and the variance of spectral estimate. This problem can be addressed by the multitaper technique [Thomson, 1982; Park et al, 1987]. 4.2.3 Mult i taper spectral analysis The multitaper method was first presented by Thomson [1982] and was applied to seis-mograms by Park et al. [1987]. The multitaper method has the lowest variance and minimum bias compared with several single-taper methods [Park et ai, 1987]. Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 91 The multitaper method multiplies the data by not one, but several leakage-resistant tapers. The statistical information discarded by the first taper is partially recovered by the second taper, the information discarded by the first two tapers is partially retrieved by the third taper, and so on [Thomson, 1982]. The multitaper spectrum, S(f), is formed as a weighted sum of several eigenspectra, , which are separately estimated from each tapered data [Park et al., 1987], N - l Yk(f) = E «n XN(n) e~i2"'n (4.70) n=0 K - l Yk(f) S(f) = £ l ^ p - (4.71) fc=0 A f e where Xj^(n) is the original data record and #(f) = E « n E ° n = 1 (4-72) n=0 n=0 I-w l#(/)f df h(N,W) = ' ' (4.73) 11% \Ak(f)\ df and k index the order of eigenspectra, a£ are the eigentapers and W is half of the spectral window width. Therefore the multitaper spectrum is already a smooth estimate; it has less variance than respective eigenspectra which have been designed to reduce bias(leakage), and it is also a consistent estimator. For the estimation of eigenspectra Vit(/), equation (4.73) can be interpreted as the fraction of energy that derives from the frequency interval | / — f'\ < W while 1 — A* is spectral energy that leaks in from outside the band. Because each taper samples the record in a different manner, the net leakage can be reduced by the offsetting effects of the different eigenspectra leakages. Only a few low order tapers can be used, as the higher-order tapers allow an un-acceptable level of spectral leakage. Also the summation of too many eigenspectra will Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 92 decrease the resolution. The preliminary results show that we should use less than 5 tapers whose A* are greater than 0.95. The arithmetic average results of eq.(4.71) can be improved by an adaptive multitaper spectral estimate which uses frequency-dependent weights, dk(f) [Thomson, 1982], * < » " A t S ^ E U U ) } ' (4-74> and produces the optimal weighted spectrum, S(f), « m EL-p1 4 ( / ) W ) | 2 where S(f) is the true value of the spectrum and Bk(f) is the spectral energy at fre-quency / that leaks in from outside the frequency band [f — W,f + W]. E{Bk(f)} can be approximated by a 2(l — Afc), where, o\ = E^o ^Ar(rc)2 is total energy of the original data. When the spectrum has a steep slope, the higher-order eigenspectra are down-weighted and the adaptive spectral estimates tend towards the least biased eigenspectral estimates [Park et al., 1987]. 4.2.3.A Evaluation: variance and bias The prolate spheroidal taper sequences are calculated with inputs of N and P, where P = NW is the time-bandwidth product of each taper. 2W = 2P/N is the frequency-bandwidth in which spectral energy concentrates. So, the spectral window of each taper is 2P x wide, i.e., it has 2P Rayleigh frequencies. In the preceding paragraphs we mentioned that variance and leakage are functions of spectral window width and the summation of different eigenspectra. Here we investigate Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 93 the role of P and the number of tapers, K, in the estimation procedure in terms of variance, leakage and total spectral energy. They are respectively defined as Var(S(f)) = j$ T,k=~o A\> multitaper [Park et al, 1987] (4.76) •7^^-—f, adaptive multitaper [Park et al., 1987] EL = I - K ^ j_ [Park et al., 1987], (4.77) E n = SkAn, ( 4.7 8 ) e JV-1 W here a\ = £ XN(n). (4.79) n=0 We conduct the investigation in a random field that is generated by the spectral method introduced in Section 4.2.3.B. We choose the variance of the random field, a\&K = 0.25, the correlation length TX = ry = 4Axw(wellbore diameter), and field size = 256x256 blocks(Axw x Axw). We estimate the spectrum with both the multitaper and adaptive multitaper methods. Figure 4.32 shows the variance, leakage and the total spectral energy with respect to P and K. From top to bottom, they are the variance Var(S(f)) ja\, the logarithm of leakage E L , and dimensionless total energy EN of the multitaper and adaptive multitaper methods. Note that by the definition in eq.(4.77), the two methods produce the same amount of leakage. The good similarity between A l and B l shows that the effects of P and K on variances of the multitaper and adaptive multitaper are almost the same. A2B2 is the logarithm leakage and it clearly shows that leakage decreases rapidly as P becomes larger. This is because more energy concentrates in [/ — + It also shows that leakage increases as K tends larger. It is because the higher order tapers are less resistant to leakages by their side lobes. Especially if P is greater than 4, this Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 94 A l A3 Multitaper method Adaptive multitaper method dimensionless variance: Var(S(f)) / o>4 1YT dimensionless variance: Var(S(f)) / o«4 Number of tapers dimensionless total energy: E N dimensionless total energy: E N Number of tapers Number of topers Figure 4.32: Variance, leakage and total spectral energy with respect to spectral window width, P = NW, and the number of prolate tapers, K, of the multitaper (left column) and adaptive multitaper methods (right column).A2B2 is the leakage for both methods. Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 95 phenomenon is much more pronounced. In A3, B3, first, we can see that Parseval's theorem is not satisfied. This is because the multitaper formulas were not required to satisfy Parseval's theorem. Secondly, we can see that the adaptive multitaper method overestimates the total energy <72, because of its attempt to compensate the energy lost to leakage by boosting the coefficients of the higher-order eigenspectra in the weighted sum, eq.(4.71). A3 tells us that the multitaper method will overestimate EN as P and K increase due to the high frequency energy gathered by the side lobes of the higher order tapers. The multitaper method underestimates EN for small P and small K because in this case the ^ weighting is not adequate. Thirdly, comparing A3, B3 with A2B2, we can see that the more the leakage, the more overestimation of the total energy. 4.2.3.B Verification: hypothetical case We test the robustness of the multitaper spectral method in two steps. First, we use a spectral method to generate a random field with a known hypothetical spectral density function. It is given by S(fx,fy) = 2 * X x X y r (4.80) where o1 is the variance of the random field and fx, fy are frequencies in x and y direction, and A x , A y are the correlation lengths in x and y direction. The random field is generated by following steps:1 spectrum—• amplitude—• complex number—• To physical space S(f*Jv) \r\ = y/S(fxJv) r = \r\co80 + i \r\sm0 IDFT random phase 9 €[0, 27r] ^ h i s is coded by Dr. R. D . Beckie, Dept. of Geol. Sci., University of British Columbia. Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 96 The parameters are chosen to be, <r2 = 1.05, A^ = A y = 18AX1", field size 512 x 512 blocks. To avoid the symmetry of the spectrally generated field, only one quarter of the field is used. We estimate the spectrum of random field by the multitaper method, and we compare the estimated spectrum with the hypothetical spectrum. Furthermore, we filter the field by a hypothetical Gaussian filter function, G(x,y) = —T—T- e *• S , (4.81) 7T Ax Ay in Fourier space it is G(fx,fy) = e-** 8 + (4.82) where A x , Ay are the filter widths. Then we also estimate the spectra of the filtered fields. We denote the primitive random field and the filtered field as Kp, Kf. We filter R~p by the Gaussian filter with different filter widths, A x = Xy = 8, 16, 32Axu' respectively and denote them as Kj, Kj6, Kj2. We estimate the PSDs by the multitaper method with different P and K parameters and we denote the multitaper estimate with P = n and K = m as PnKm. We compare them with the periodogram and the true spectrum. The results are illustrated in Figure 4.33. Note that in Fig. 4.33, we have averaged the 2D spectra over the frequencies / x , fy at distance / r , to reduce the dimension to ID, i.e., the ID spectra are the averaged slices of the 2D spectra, S(fr) 1 2 A 2 A = M £ £ ^ fx—O fy—O where M is number of data points averaged. Furthermore, we scale the frequency by the the grid block scale Ax, i.e., / * Ax, resulting in dimensionless frequency. For simplicity, the 'dimensionless' is omitted. All frequency space plots in the rest of the thesis are graphed in dimensionless frequency. Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 97 target spectra of primitive field and filtered fields and Rayleigh frequencies o.i FREQUENCY primitive field le+03 lo+02 l e f O l le+00 le-01 le-02 le-03 le-02 FREQUENCY le-01 filtered by X = 8 le+03 le+02 lo+Ol k l&fOO t le-01 le-02 b le-03 FREQUENCY continued on next page • Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 98 Figure 4.33: Comparisons of multitaper spectra, periodogram and true spectra. The primitive field is filtered by a Gaussian function G(x,y)=^ A 6 A exp[-6(f^- + A4)], where filter width A x = \y= 8, 16, 32Ax. Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 99 Of the estimates for the primitive field, the periodogram and the multitaper spectra P2K1, and PsK3 match very well with the true spectrum S(f). However, P^QKQ deviates from S(f) at frequencies smaller than 0.1. For the filtered field Kf, the P2K\ and P5K3 still fit the true spectrum, but the P30K6 deviates further from S(f) and the periodogram deviates S(f) at / > 0.1. The PZQKQ and periodogram of Kj6,Kj2 continue to deviate from S(f). The P2KX of Kf starts jumping above S(f) at / > 0.09 while in Kf the deviation begins at / > 0.03. The primitive and the filtered fields are 256 by 256 pixels each. So the Rayleigh frequency is ^ j . The top graph in Fig. 4.33 shows all the true spectra and the Rayleigh frequencies of 2D data series with N by TV data points, where TV = 2 3,2 4,2 s, 2 6,2 7,2 8 respectively. Rayleigh frequency ^ produces the best resolution and the resolution of K^4 turns out the poorest. 4.3 Testing the Wiener filtering method We test the system identification technique in Section 4.1 by using it to invert the hypo-thetical Gaussian filter function that was used in Section 4.2.3.B in filtering the primitive random field Kp. We use the Kp and the filtered fields Kj as the system input and out-puts. The PSDs estimates are as those of Section 4.2.3.B. We use the iterative procedure of Fig. 4.30 to identify the filter function. Figure 4.34 shows the estimated filter functions in spatial and frequency domains against the hypothetical Gaussian functions. Note that the filter functions are averaged over the data points (x,y) at distance r, = jfx £ £ o(.,») i w l x=0 y-0 1 7V-1 N-l A *—> V—> (4.84) Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 100 and the transfer functions are averaged over the frequencies fx,fy at frequency fr, 1 l 1 2A 2A ^ '"VfM = W2 ^ (4-85) 2 /*=0 /„=0 where Afi, M2 are number of data points averaged. We can see that for \ x = Xy = 8,16Ax w , the Wiener filtering method almost backs out the hypothetical function exactly. The misfits for Xx = Xy = 32Axw are because the Rayleigh frequency is too coarse for the resolution of the spectrum of Kj2, see Fig. 4.33. For better estimation, we need increase field size to decrease the Rayleigh frequency. Chapter 4. SYSTEM IDENTIFICATION: SPECTRAL ANALYSIS 101 HYPOTHETICAL ESTIMATED -«— FILTER WIDTH = 8 • i i i HYPOTHETICAL -ESTIMATED -FILTER WIDTH =8 4 2 3 4 5 6 7 8 PHYSICAL SPACE 0.10 0.15 FREQUENCY HYPOTHETICAL -ESTIMATED • FILTER WIDTH = 16 HYPOTHETICAL • ESTIMATED -FILTER WIDTH = 16 5 10 PHYSICAL SPACE HYPOTHETICAL ESTIMATED -•— FILTER WIDTH - 32 1 l **"•—»• i i HYPOTHETICAL ESTIMATED -*— -FILTER WIDTH = 32 1 1 * T * ~ t l i l • 1 i i 10 13 PHYSICAL SPACE 0.00 0.01 0X2 003 0.04 0X5 FREQUENCY 0.06 0.07 0.08 Figure 4.34: Left column: comparisons of the estimated filter function with the hy-pothetical filter function G(x,y) = _ 6 ( + -V ) e x* xy ; right column: com-parisons of the estimated transfer function with the hypothetical transfer function G{fx, fy) = e ~ e *2 ( A ' + ). w h e r e filter w i d t h X x = X y = g ) 1 6 ) 32Ax respectively. Chapter 5 S Y S T E M I D E N T I F I C A T I O N : P E R T U R B A T I O N M E T H O D In this chapter, we utilize a numerical perturbation method to identify the system filter. A general introduction to perturbation methods is given in Section 5.1. The mathe-matical formulas of the numerical perturbation method are derived in Section 5.2. The perturbation method requires the inverse conductivity estimation, already discussed in Chapter 3, from the response curve measured in a perturbed conductivity field, and the accuracy of the method is dependent on the inversion procedures. Thus, in Section 5.2.1 one more inversion scheme is analyzed and compared with the scheme a and 6 of Chapter 3. The area and the amount of perturbation are very crucial for the accuracy of the per-turbation method. The area and amount of perturbation which possess the minimal bias and the highest accuracy are numerically determined in Section 5.3. The system filters with different specific storages and the appropriate perturbation schemes are discussed in Section 5.4. 5.1 Introduction Perturbation methods are widely used to identify the statistical relationship between the fluctuations of target variable, <f>(x), with the variations of its input parameter, p(x). Both <f>(x) and p(x) are often seen as random processes that can be expressed as an expectation or mean, plus a perturbation, i.e., <f>(x) = E[<f>(x)] + e^ , p(x) = E\p(x)] -f ep, where and ep denote the perturbation and they are required to be relatively small compared to their mean values for the convergence of the perturbation. 102 Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 103 In stochastic subsurface hydrology, perturbation approaches are extensively employed to derive the statistic dependence of hydraulic head variance on the variance of logarithm conductivity and the relationship between effective conductivity and the small scale het-erogeneity [Gelhar, 1994]. The assumption of a constant mean hydraulic gradient, i.e., uniform mean flow in space, greatly facilitates the perturbation analyses [Gelhar, 1994]. However, this condition may not be met in reality, for instance, in the nonuniform radial flow to a pumping well. In the nonuniform flow study, Oliver [1990] characterized the relationship between drawdown and permeability as a Frechet kernel function of the radial flow in pumping tests, Si = J Gi(x) (1 - dx; KD(x) = (5.86) where Gi(x) is the kernel and (1 — K^x^) is the presumed permeability distribution which is radially symmetric but nonuniform, and K is the average permeability. He solved the kernel by the perturbation SD = PD0+^PD1+e2PD2+O(ts) (5.87) where the perturbed drawdown SD is approximated by the sum of pressure terms of order smaller than 3, and e is a small perturbation. This kernel function was used by Oliver [1992] to derive the permeability distribution in eq.(5.86). Oliver [1993] gener-alized the (1 — K*(x)) permeability pattern to nonradially symmetric porous media and analytically solved the associated Frechet kernel by a perturbation method similar to eq.(5.87). Another method similar to the perturbation idea is the sensitivity analysis, which can be regarded as a special case of the perturbation approach in which <f>(x) and p(x) usually are deterministic, i.e., = ep = 0. Suppose the functional dependence of <j>(x) Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 104 on p(x) is in the form <f>(x) = g(p(x)), (5.88) then for variation of p(x) about some value, say its mean pp, denoted as 8P = p(x) — pp, <l>(x) can be approximated by a series in terms of pp and 6P, for instance, the Taylor expansion is The d9QpPffl term in above equation is the sensitivity of p(x) to <f)(x), the associated analysis is referred as sensitivity analysis. Yeh [1986] reviewed several commonly used sensitivity methods. Wang [1992] used a sensitivity method to identify the principle contamination sources to the shallow aquifer underlying a chemical-industry city where various polluting factors had been detected. The groundwater quality management model of Dong [1994] is also an application of sensitivity analysis, i.e., the objective function is calculated from the response matrix that results from the variation of the flux and concentration under planning. 5.2 Numerical perturbation formula The pros-and-cons of analytical perturbation versus numerical ones are broadly discussed by Gelhar [1994]. The chief disadvantage of analytical methods is their dependence on oversimplified conditions which may be far from reality. Numerical methods can account for practical situations in much greater detail. In the meantime, the generic meaning of numerical results could be as good as the analytical ones if appropriate numerical methods are used. Here the task is to use a perturbation method to deconvolve the system filter from Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 105 eq.(3.29), which is Ym(x) = r G(x-x')Yc(x') dx' + N(x) (5.90) J — CO where Ym(x) is inverted from the response curve H(t) measured in core-scale field K:(#)-N(x) accounts for model misfits mainly due to discretization errors. We shall use the high accuracy 2D FDM model for the simulation, 1 thus the noise term becomes negligible and the above equation simplifies to /CO G(x - x') Yc(x') dx'. (5.91) -co One may attempt to apply analytical perturbation method and directly derive H(t) for a given Yc(x), 2 and then use the H(t) solution of Cooper et al. [1967] to obtain Ym, at last the resulting functional dependence between Ym and Yc(x) gives the filter G(x), i.e., Yc(x) — • H(t) — • Y m ( x ) 4 I G(x) But the inverse step H(t) —• Ym(x) is very complex to derive analytically, 3 as can be seen from the H(t) equation, 8 H0 a f°° du H{t) = — / e ° — , 5.92 V IT2 Jo u • Au v r2 10Ym • t a = — Ss, P = rl Au = [u J0(u) -2 a Jx(w)]2 + [u Y0(u) -2 a Y^u)]2. The complexity of above equation dictates the impossibility of analytically obtaining the functional relationship Ym ~ H(t) without some simplifications, that may in turn 1 W e have checked its accuracy in Chapter 2. 2Perform perturbations on Yc(x), like that of Oliver [1993]. 3 W e have discussed the numerical inversion of H{t) —• Ym(x) in Chapter 3. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 106 cause intolerable errors so as to plague the analytical attempts. Contrarily, the steps Yc(x) —• H{t) and H(i) —• Ym(x) are convenient to be performed numerically. With the numerical perturbation, Ym(x) can be readily obtained from the inversion estimation of RC (response curve), and Yc(c) is known, and what we are after is the spatial filter G(x). Let us first discuss the simple case of homogeneous Yc(x), in which Ym(x) = Yc(x), then we get /oo G(x-x')dx' = 1 (5.93) -oo it means the filter volume equals to 1 for slug test at location x. But we still don't know the filter amplitude distribution around x, and for that purpose we must lift the integral /f^,, which is the deconvolution task of the numerical perturbation, discussed in the following. Figure 5.35: The radial mesh for the simulation of perturbation experiments. We perturb Yc(x) by an amount of 8yc in an area of Clp, Fig.5.35, 4 the change in Ym(x) can be written as 4 Here the numeric perturbation scheme is slightly different from the traditional analytical method. Analytical perturbation is usually performed in the whole parameter domain by a given statistic de-scription of the perturbation, e.g., uep, c\ . Here we numerically perturb, only once a time, in a small area fip in the parameter space. This can be achieved easily by a discrete numerical method, while in contrast, such a perturbation is onerous to be solved analytically. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 107 6ym{x) = (Dym,8Yc) + 0[6Yc(x) (5.94) where 8Ym is the change in measured conductivity from the perturbation 8Yc, and (DYm, 8Yc) is the Frechet derivative G(x - x') 8Yc(xl) dx', -oo (5.95) where 8YC is >yc(*) = < 0, x g ftp, constant, x € ftp. (5.96) 0 [^ yc(x)] represents the high order terms. Note that eq.(5.95) is actually is the notation for (Dym,8Yc) = f°° G(x-x') \Yc(x) + 8yc{xl)} dx' J—oo L J /oo G(x - x') Yc{x) dx'. -oo (5.97) substituting eq.(5.91) to eq.(5.97), we have Ym(x) + SYm{x) = G(x-x') [Yc + 8Yc{xl)] dx' + O [^ J . (5.98) This is the initial perturbation equation. Clearly, the filter function is dependent on the flow pattern of slug test 5 and the perturbation will alter the local flow pattern. Thus the filter function is also being perturbed. The perturbed filter function, Gp(x), can be written as Gp(x) = G(x) + 8G{X) (5.99) 'Because Ym is inverted from RC and is dependent on the flow pattern. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 108 where &G(X) is the filter variation invoked by the b~Yc(x) perturbation. Thus the high order term in eq.(5.98) equals = r tG(x-x>) [Yc(x) + b-Yc(xl)] dx' (5.100) J —oo If we assume 6G(X) = 0, i.e., O t>Yc(x')} = 0 so that eq.(5.94) reduces to /CO G{x-x')6Yc(x>)dx' (5.101) -co &YC(X) is constant over fip, denoted as b~Yc\np z e r o everywhere else, eq.(5.96), thus we have = / G(x - x') 6Yc dx' UP = *y«lnp / G(x-x') dx' (5.102) which is / G(x-x') dx' = (5.103) If we further assume that G(x) in fl p to be approximately uniform, denoted as G\np, then (5.104) G\n « J- s-^l P ° f sYc\n„ This is the final perturbation formula, meaning that: the filter amplitude in ilp = ( the change of the conductivity measured at wellbore located at x in a homogeneous aquifer where the conductivity in fip is perturbed by £y c | n p ) -f-[(the perturbation magnitude 6yc|n p ) x (the perturbation area Qp)]. In the above derivation, there are two operations on the filter integral, 1) Move the intergral from [—oo,+oo] to fip, eq.(5.91) — eq.(5.103) G(x - x>) dx' —» /flp G(x - x') dx' 2) Assume G is constant in f l p , thus taking it out from under the integral. fQp G(x - x') dx' — G\ap f n p dx = G\nnp Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 109 The assumptions made for each operation are: 1) Filter function does not change under perturbation, i.e., SG(X) = 0. 2) G(x) is approximately uniform in fl p . These assumptions can be asymptotically approached by proper choice of 6yc and ftp, which is crucial for the perturbation equation (5.104) to be valid, and which is very complex to be determined analytically. We use numerical techniques. In Section 5.2.1 we compare different inverse estimation procedures of perturbed con-ductivity, and in Section 5.3 we employ numerical techniques to determine the area and the magnitude of the numerical perturbation. 5.2.1 T h e inverse estimation procedures of perturbed conductivity We define the perturbed logarithm conductivity, Yp, as Yp = Yc + 6Yc (5.105) = \g(Kc) + lg(«) where 6yc is the logarithmic perturbation. Thus the perturbed conductivity, Kp, is Kp = aKc. (5.106) We define the dimensionless perturbation amount as 4 = \I^IA = |! _ fl|| ( 5 . 1 0 7) which represents the extent of perturbation. If we perturb dk, once a time, at the blocks along the same line shown in Fig.5.35, then there may be three inversion methods for the RC measured for each perturbation: a) Estimate only IgK with lgS^ equal to true value; b) Concurrently estimate lg.ft' and lg S„; c) Estimate lg K at each time step with lg Ss equal to true value. We have discussed method a and b for the RC measured in heterogeneous conductivity fields. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 110 Now we compare the estimates by method c with those of a and b. 6 Here forth, we use superscripts °, 6 , and c to denote the estimates from these methods. The results are shown in Fig. 5.36 where the unperturbed lg K and lg S3 equal to 2 and —1 respectively, and the perturbation dk is 0.5. 7 With method a, Fig. 5.36A shows \gKa from each perturbation versus the corre-sponding dimensionless perturbation locations rD = r/rw. Because dk equal to 0.5 is a negative perturbation, the measured lg K is less than the unperturbed lg K. Away from wellbore, the perturbation area ftp increases(the perturbation strength, approximated by ftp X dk also increases), while the influence of the inhomogeneity on RC decreases, the net effect results in the curve shape shown in Fig. 5.36A. With method b, graph B and C of Fig. 5.36 demonstrate \gKb ~ rD and lgSb ~ rD. We can see that \gKb is biased because of the deviation of lgS^ from the true value, and wherever lg Kb is underestimated, lg S^ is overestimated, and vice versa. Thus we do not choose this method, as we discussed in Chapter 3. We also discussed in Chapter 3 that method a may possess bias because the perturbed wellbore response curve can not be entirely matched by that of a homogeneous model, but this error is far smaller than that from method 6, and also it can be minimized by an optimal choice of a minimum nontrivial 6yc and a proper Clp. 8 With method c, Fig. 5.36D displays lg K° ~ rD where lg K° at rD is the arithmetical average of the lg at rD inverted at all time steps, i.e., the geometric average of Kc.9 6Here we need to discuss these inversion methods because they result in different lg K estimates, hence, the different filter functions associated with these estimates. We did not use method c in Chapter 3, see the arguments therein, but here for the inversion of the RC of the theoretical perturbed conductivity fields, the estimates from method c can be readily obtained, and they can be used to calculate the filter function at each time step or for all time at certain spatial location. 7We proved that in Chapter 3 that the perturbation is independent of the unperturbed \gK, thus we can choose any value without losing generality. The other lgS, cases will be discussed in Section 5.4. 8Another reason to reject method b is that in the perturbation case the 6Ym and the 6yc usually are small compared to the unperturbed Ye, and the small changes of Ym due to 6yc may be severely damped by the interactivity of Y^ ~ lg Sj in the inversion, which will be shown in Section 5.3. 9The choice of such an average seems more arbitrary than with some theoretical proof. The possible Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 111 The comparison of above lg K estimates are shown in Fig. 5.36E. We see that lg Ka lg Kb intersect at the point P2 of graph B, where rD = 1.8. This is the only point at the influence zone where the b estimation is unbiased. . ' The lgKa intersects the \gKc at two locations, where rD = 1.3, 5.1. If we use lgKa as the standard and evaluate the deviation from it as under- or overestimation, then \gKc is overestimated at [0, 1.3) and (5.1, 50], underestimated at [1.3, 5.1]. Note that method b and c overestimate the influence radius Re,10 where R^ = 30, R% = 50, whereas Rae = 19.5. - -The lg Kh also intersects the lg Kcat two points. Interestingly, the most overestimated point of method c overlaps the most underestimated point of method b, where rD equals to 1.1. This radius might be interpreted as a "numeric-estimation wellbore skin" radius, in the sense that the estimated conductivity at such radius has abnormal high or low values. Note that if a positive perturbation is used, the lg Kb will be most overestimated at rD = 1.1, as we illustrated in Chapter 3. The above IgK estimates give rise to different filter functions that are shown in Fig. 5.36F. Their intersections are the same as that in graph E. If we do not average the estimates from method c, we can appreciate the transient filtering property from two aspects: Firstly, illustrated in Fig. 5.37A, at each point in the influence zone, the filter function changes with time. At locations closer to wellbore, the filtering has a higher power and happens earlier. The filter function at all locations decay off to zero until the end of total time of which /3 = Kt/r^ = 50. Secondly, shown in Fig. 5.37B, the filter function in dimensionless space rD G [1, Re] varies with time. At earlier time most of the filtering power concentrates around wellbore and theoretical consideration is that the temporal average is co-phased with spatial average, which is reported in various literature that a geometric average is appropriate. 1 0The concept of influence radius was discussed in Chapter 2, here the influence radius is scaled to be dimensionless by rw. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 112 log10(K")~rD zoom 2.0000 "1 1.9999 1 g 2 1.9998 1.9997 1.9996 }• perturbed log l0(K) unperturbed loglO(K) 10 P e r t u r b a t i o n l o c a t i o n 100 logio(K^-rD »oglo(Ssc) - r D 2.0002 2.0000 1.9998 2 19996 1.9994 -0.9985 -0.9990 | -0.9995 BO •a J -1.0000 e. o J"-1.0005 -1.0010 -1.0015 perturbed log l0(K) unperturbed log l0(K) overestimate logio(K) between PI and P2 underestimate logio(K) between P2 and P3 PI P2 P3 B 10 Perturbation location 100 perturbed logl0(Ss) unperturbed loglO(Ss) P3 underestimate loglo(Ss) between PI and P2 overestimate logio(Ss) between P2 and P3 le+Ol Perturbation location le+02 continued on next page Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 113 logiotfT) ~ r D 2.0001 2 . 0 0 0 0 I? "5 1.9999 S 1.9998 1.9997 1.9996 per turbed l o g l O ( K ) unper turbed l o g l O ( K ) D 10 Per tu rba t ion l o c a t i o n 1 0 0 comparison of logio(K')~r D l o g i o ( K V r D log 1 0(K<)~rD 2 . 0 0 0 2 2 . 0 0 0 0 1.9998 1 1.9996 1.9994 i n v e r t l o g 1 0 ( K ) o n l y i n v e r t l o g l 0 ( K ) w i t h l o g l 0 (Ss ) inve r t l o g ( K ) at each t i m e step unper turbed l o g l u ( K ) Pe r tu rba t ion l o c a t i o n comparison of filter functions G ' ~ r D G b ~ r D G c ~ r D 0 . 3 0 0 . 2 0 0 .10 fi CL. 0 .00 - 0 . 1 0 - 0 . 2 0 - 0 . 3 0 i n v e r t l o g l 0 ( K ) o n l y inve r t l o g l O Q C ) w i t h l o g l 0 ( S s ) inve r t l o g ( K ) at each time step F Per tu rba t ion l o c a t i o n 10 Figure 5.36: The lg K inverted by different methods and the associated filter functions. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 114 Figure 5.37: Spatial filter is a function of time. A shows the filter function at locations Pi, P2, P3 at time [0,3], 8 = 50. B shows the filter function at time ti, t2, t3 in dimensionless space rD, rD G (1, Re]. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 115 decreases rapidly, while at later time the averaging power spreads to a larger volume and its amplitude drops off much more slowly. For the same reason discussed in Chapter 3, we shall invert the lg K with the whole response curve by method a in following discussions. 5.3 T h e area and the amount of perturbation Analytical methods allow one explicitly determine limits o f perturbation through evalu-ating the approximations used in the perturbation and the final analytical expressions. For the numerical perturbation, however, we don't have explicit results, thus we have to decide the area and the amount of perturbation by a "try-evaluate-try" approach that is guided by the proper interpretation of the intermediate numerical results with respect to the physical behavior of the slug test system. From eq.(5.104) we know that, physically, ftp can not be too large otherwise G(x) is not close to be uniform in ftp and such a perturbation may strongly alter the flow pattern so that it invalidates the assumption of eq.(5.101); Also, ftp can not be too small otherwise the perturbation becomes trivial; mathematically, if ftp —> 0, then G(x) —• +oo, it causes overestimation; if ftp —• +oo, then G(x) —> 0, it introduces underestimation. We select ftp by two steps. Firstly, we adjust the block sizes of the 2D radial FDM model of a slug test in a homogeneous aquifer until the numerical RC matches the RC of analytical solution. This means the space discretization is suitable to resolve the hydraulic gradient in the influence zone. The same criteria is used to choose ftp, i.e., a suitable ftp is the one in which hydraulic head is approximately uniform. So secondly, we simply select blocks along one line, opposite to which is a nonflow boundary. This nonflow boundary is for the angular variable 0, and it is physically valid because the Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 116 perturbation invokes symmetric flow pattern that warrants no crossing flux. The perturbation dk can not be too large, otherwise it alters flow pattern drastically, e.g., one limiting case, if Ss = 1, and dk = 100, it turns the perturbed block into an open well. Preliminary tests suggest it is better to choose dk < 1, i.e., 0 < a < 2. On the other hand, dk can not be too small in which case the changes of the perturbed wellbore response curve might be close in magnitude to the tolerance of the numerical equation solver, hence the calculated filter function might be instable. Furthermore, if dk too small, the influence of perturbation at blocks farther away from wellbore would not affect wellbore response, thus resulting in an underestimation of the volume and averaging power of the filter function. We find the optimal perturbation amount through a number of tests of different dk values. Two criteria should be met for an appropriate perturbation: 1) The integral of the filter approximately equals to 1, 1 1 2) The skewness —* 0, where K* = (Ym of negative perturbation) + (Ym of positive perturbation). Under this condition, the same filer function G\ap « gyl^ should be obtained no matter what kind of perturbation is performed. Fig. 5.38 shows the results, estimated with method a and b, of dk = 0.5, 0.1, 0.05, 0.01, where 0 < a < 1, unperturbed lg K = 2, and unperturbed lg Ss = — 1. 1 2 Fig. 5.38A and Fig. 5.38B show lgA'" ~ rD and the associated filter functions Ga ~ rD. We see that the filter functions converge to the one of dk = 0.01, and we shall show in the next section that the filter integral approximately equals to 1. Thus we choose the optimal dk 1 1 Because the perturbation is usually only a small fraction of the unperturbed homogeneous conduc-tivity, the the filtering should be very close the homogeneous case. The filtering of a unit filter volume is called mean-preserving filtering. 1 2 Here we use negative perturbations. In Chapter 3 we illustrated that the skewness K™ —+ 0 if a > 0.9, i.e., dk < 0.1, Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 117 2.0001 Figure 5.38: Perturbation with different <4- The IgK is estimated with lg,Ss equal to true value. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 118 measured logio(K) of different perturbation Per turba t ion l o c a t i o n 1 0 0 -0 .9990 measured log10(Ss) of different perturbation dk filter function of different perturbation - 1 . 0 0 1 0 0 .05 0 . 0 0 -0 .05 -0 .10 -0 .15 -0 .20 -0 .25 h Per turba t ion l o c a t i o n • 1 d k = 0 . 5 0 -k d k = 0 . 1 0 dk=0 .05 •V dk=O.01 1 // f I c i i i • —I i i 1 • 1 • - 1 Perturbation location 10 Figure 5.39: Perturbation with different dk. The \gK is concurrently estimated with lgSa. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 119 as the one to which the filter functions, estimated with a range of dk values, converges to 1. 1 3 Here the optimal perturbation amount dk is 0.01. In contrast, with method b the \gKb ~ rD and the associated filter functions Gb ~ rD appear inconvergent, shown in Fig. 5.39. This inconvergency is caused by the interactive property of lgKb and IgS^ estimation that we analyzed in Section 5.2.1. This further justifies our judgment that method b is not appropriate as the inversion procedure for the perturbation. Clearly, for different specific storages, the perturbation affects wellbore response in a different manner. Smaller Sa has stronger influence at farther locations while larger Sa is more sensitive around wellbore. Thus the selection of dk should be S -^dependent. So it is no wonder that a suitable dk for a small Sa may be too large for the perturbation with a large Sa. The Sa in Fig. 5.38 is 0.1. For Sa < 0.1 we may need a dk larger than 0.01. This is the topic of the next section. 5.4 Fi l ter functions of different specific storages We use several statistical measures to characterize the filter functions of different specific storages. The volume of filter function, Va(rD) is defined as Vo(rD) = J 2irG(rD)rDdrD, (5.108) where rD is dimensionless distance, and £ is a infinitely small positive number, i.e., ^ —>• 0 +. Because the 2D radial FDM we used is mesh-centered, the first node is on the inner boundary where rD equals to 1, so that G(l) is zero and G(l + £) represents filtering power at the node closest to wellbore in the aquifer. Re is the dimensionless influence 1 3In the lg Ka ~ rD plots, we can observe that the curves of different dk values converge to a curve of some dk value, and filter integral associated with the dk should approximately equal to 1. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 120 radius defined as d[VG(rD)} 0, (5.109) d(rD) which means the filter volume stabilizes at and beyond Re. This may be a more rigorous definition than that by Guyonnet et al. [1993]14 because: 1) the ^Q I^^ ^ explicitly incorporates the concept of gradient which manifests the rate of increase of the filter volume; 2) the G(rD) in the expression of [Vcr(rr>)] explicitly tells the contribution of the conductivity at rD to the conductivity inverted from the wellbore response curve. The dimensionless equivalent filter width is defined as K = ( M l . ) In the next Sections 5.4.1 and 5.4.2, we shall compare the filter functions that are estimated with two perturbation schemes: 1) Fixed perturbation magnitude, and 2) Ss-dependent perturbation. 5 . 4 . 1 Fixed perturbation magnitude We fix the perturbation magnitude, <4 = 0.5, for Ss = 10~n, n = 1, 2, . . . , 10. Con-sequently, from the discussions in the previous section we expect filter functions to be slightly overestimated for large Ss and slightly underestimated for small storage. To distinguish the different filter shapes, the filter functions are plotted in logarithm coordinates, lg^r^)] ~ lg , in Fig. 5.40A. It displays that, as Ss decreases, the filter support area increases, and the amplitude of the filter away from the wellbore decreases less rapidly. Fig. 5.40B is a plot of Va{rD) ~ l g r D - ^ s r i o w s 1) VG(rD) ceases increasing as rD approaches Re, 2) With larger Ss, Va(rD) grows faster while the influence area be-comes smaller, and 3) Va{Re) approaches a unity volume for Ss decreasing from 10 _ 1 14Their definition is "the maximum distance traveled by a perturbation ^ = 1%, 5%, 10%". Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 121 to 10 - 8. The volume overestimation(VG(i2e) > 1) is what we have predicted, e.g., in Fig. 5.38Al(where Sa = 10_ 1), the filter function of dk = 0.5 is overestimated. The Va(Re) of Ss = 10 - 9, I O - 1 0 is less than 1. This indicates the perturbation of dk = 0.5 is too small so that the filter functions are underestimated. Fig. 5.41A demonstrates the linear relationship between Ae and lg(^ 5=)- Note that since the overestimation or underestimation of Vo(Re) and G(l+£) occurs concurrently,15 thus = Ae is a more reliable representation of the filtering property than the G(rD) and Va(rD) of Fig. 5.40 are. Fig. 5.41B is a plot of lg(-Re) versus lg(^ y=). The Re is taken from Fig 5.40B. It shows a good linear relationship. 5.4.2 S*-dependent perturbation Following the same reasoning of Fig. 5.38A1, we perform the perturbation tests of dif-ferent dk for Sa = 10 _ n, n = 1, 2, 10. Some of the test results, G(rD) ~ rD and ~ ^ § ( ^ 7 ) 5 a r e graphically shown in Fig. 5.42 - 5.45 whose Sa equals to 10 - 1, 10~3, IO - 5 , 10 - 7 respectively. The curve of dk = 0.5 can be used as a reference for comparison across these graphs. In the graph A of Figures 5.42 - 5.45, the optimal dk values are found to be 0.01, 0.05, 0.1, 0.3, respectively. The smaller the Sa, the larger the optimal dk. This tells us that aquifers of high specific storage are more sensitive to the conductivity variation. In Fig. 5.42, 5.43, we can see that, for the perturbations smaller than the optimal dk, the filter functions become instable because the hydraulic head changes from perturbation are likely of the same magnitude as random truncation errors of the matrix solver. In the graph B of Fig. 5.42 - 5.45, we can see that the total averaging power of the filter function calculated with the optimal dk approaches unity. We also can see that for 15Because Va(Re) is an spatial integral of G(l + £)• Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 122 the filter functions calculated with dk greater than the optimal value, the greater is the dk, the more overestimation is the Va(Re)', Note that they have identical influence radius; This conforms to the reality that the influence radius is independent of conductivity and it is only a function of specific storage. Note that, in graph B of Fig. 5.42 - 5.45, we have approximated the integral of VG(rD) defined in eq.(5.108) by two steps. We first interpolate the estimated filter functions into 106 data points evenly log-spaced; Secondly we sum up these discrete points. The associated interpolation error could contribute both positively and negatively to the averaging power Vb(rD), depending on the shape of the filtering functions, but in general they should be insignificant. Fig. 5.46A shows the Va(rD) ~ lg(^-) calculated with the optimal dk for Sa = 10 - n, n = 1, 2, . . . , 10. Fig. 5.46B shows the comparison of the Ae of Fig. 5.46A with the Ae of Fig. 5.40B estimated with dk = 0.5. This graph shows 1) Ae ~ s^(^ =) obeys a linear relationship, 2) the filter width of the filter functions calculated by dk = 0.5 are more underestimated for decreasing Ss. This can be explained with the Vo(rD) ~ lg(^=) pattern in Fig. 5.40 and by the Ae definition in eq.(5.110). Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 123 le+OO le-02r le-04 I E le-06t le-08f le-10 le-13-le-14 le-16* le+OO le+01 le+02 le+03 le+04 Dimensionless distance Ss=le-10 le+05 le+06 1.6 1.4 1.2 1 1 s 3 0.8 o u •3 0.6 > 0.4 0.2 •Ss= Influence radii 1.95el 6el 1.95e2 6e2 1.95e3 6e3 1.95e4 6e4 1.95e5 6e5 Ss=le-10 le+OO le+01 le+02 le+03 le+04 Dimensionless distance le+05 B le+06 le+07 Figure 5.40: The spatial filter for different specific storages, Sa = 10 _ n, n = 1,2,..., 10. For all perturbations, <4 = 0.5. A shows the filter function versus dimensionless distance in log space. B is the volume of filter function integrated over dimensionless distance r , Re VG(TD)= f 2irG(rD)rDdrD. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 124 50 40 30 20 10 le+01 le+02 le+03 le+04 Specific storage! l/sqrt(Ss)] le+05 le+06 le+01 le+01 le+02 le+03 le+04 Specific storage[l/sqrt(Ss)] le+05 Figure 5.41: A is the equivalent filter width defined as Ae = • B is the influence radius Re defined as —* 0. 9 M Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 125 Figure 5.42: Sa = IO"1. A shows that G(rD) ~ rD converges to the dk = 0.01 curve. If dk is too small(c?fc < 0.01), solutions become instable. B is the corresponding Vo(rD) ~ Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 126 Figure 5.43: Ss = 10 3 . A shows that G(rD) ~ rD converges to the dk = 0.05 curve. If dk is too small < 0.05), solutions become instable. B is the corresponding VG(rD) ~ Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD Figure 5.44: Sa = 10 5 . A shows that G(rD) ~ rD converges to dk = 0.1 curve. B is corresponding VG(rD) ~ lg(^;) curves. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 128 1.2 le+OO le+01 le+02 le+03 le+04 le+05 Dimensionless distance Figure 5.45: Sa = 10 7 . A shows that G(rD) ~ rD converges to dk = 0.3 curve. B is the corresponding VG{TD) ~ lg(^-) curves. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 129 continued on next page Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 130 le+OO le+OO le-rOl le+02 le+03 le+04 le+05 le-f06 Dimensionless distance Figure 5.46: A shows the Vo(rD) ~ lgr^ curves calculated with the optimal dk for Ss = 10~n, n = 1,2,..., 10. B is a comparison of the Ae ~ lg(^j) of A with those estimated with dk = 0.5 for all S, values. C is the corresponding lg [G^)] ~ lgr D . Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 131 5.4.3 Parametric representations In this section, we present the parametric representations of the filter functions, equivalent filter widths and influence radius, determined from the~6ptimal dk for different specific storages in Section 5.4.2. The filter functions lg[G(rD)] ~ lg*^, shown in Fig. 5.46C exhibit straight-line patterns in the intermediate portions which constitute the major parts of whole curves. Here we only fit the straight-line portions by a parametric representation, G{rD) = G(l+0r4, ^ 0 , (5.111) VD) where the G(l + £) is the filter amplitude at the wellbore edge, and p is an exponent. Table 5.4: The G(l + £), P of the filter expression G(rD) = G(l + £) -r-^. Sa IO"1 IO"2 10"3 10"4 IO"5 10"6 IO"7 10"8 IO"9 l f J-io G(l+£) x IO"2 18.338 7.908 4.882 3.472 2.840 2.412 2.165 1.950 1.885 1.753 P 3.0 2.3 2.1 2.05 2.0 2.0 2.0 2.0 2.0 2.0 Table 5.4 shows 1) G(l + £) decreases with decreasing Sa, and 2) G(rD) ~ rD approxi-mately approaches a . \ a decreasing law for the Ss smaller than 10~4. Now we fit both {RD) Gil + 0 ~ T^T a n d p ~ v ^ T ' s n o w n m F l S - 5 - 4 7 - T n e fitted c u r v e s a r e G(l-rC) = 0 .25 0 2 5 + 0.2S°S-75 + 0.01753 (5.112) p = 3 . 0 ^ + 2.0 (5.113) Thus eq. 5.111 becomes G(rD) = (0.2Sr + 0 .2S 0 7 5 + 0-01753) ^ . (5.114) \rD) The above equation describes the direct relationship between G(rD) and Ss. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 132 0.25 czt 0.20 x +• 0.15 ff u f 1 I 0.10 I S 0.05 Fitted G(l+xi) versus Ss G(l+xi) versus Ss 0.00 1 le+01 le+02 le+03 le+04 le+05 Specific storage l/sqrt(Ss) le+06 le+07 3.00 2.80 p. 2.60 1 8-2.40 2.20 2.00 Fitted p versus Ss p versus Ss B 1 le+01 le+02 le+03 le+04 le+05 Specific storage l/sqrt(Ss) le+06 le+07 Figure 5.47: A: Fitting G(l + £) ~ by G(l +() = 0.2S S ° 2 5 + 0.2S5075 + G.01753; B: Fitting p ~ ^ by p = 3.0y/Ss + 2.0. Chapter 5. SYSTEM IDENTIFICATION: PERTURBATION METHOD 133 The approximate straight-line portion of A e ~ lg(^-), shown in Fig. 5 .46B, can be parameterized as Ae = A° - [ l g ( ^ ) ] ° + (5.115) A 0 = 5.47, [ l g ( ^ ) ] ° = 0.5 (5.116) The approximate straight-line portion of R^ ~ lg(^j)> shown in Fig. 5.41, can be parameterized as lgi?e = (lgi? e)° - P g ( ^ r ) ] ° + U 5 ( ^ ) ( 5 . H 7 ) R e = (5.118) ySa (lgi? e)° = 1.29, [ l g ( - ^ ) ] ° = 0.5 (5.119) Chapter 6 F I L T E R F U N C T I O N In this Chapter we investigate the filter functions in heterogeneous core-scale fields of various spatial variabilities and specific storages, utilizing the Wiener filtering method-ology. To explore the effects of Ss on the spatial filtering, in Section 6.1 we perform slug tests in isotropic Kcore fields of 64 x 64, 128 x 128 and 256 x 256(Ax"' x Axw) with S3 — 10 - n , n = 1, 2, . . . , 5. In Section 6.2, we present the effects of the spatial variabil-ity of conductivity on the spatial filtering. We consider the spatial variability of variance, correlation length, variogram model and distribution pattern. The investigations in this chapter are not designed to accommodate every cases met in practice, but intended to disclose the general roles of Sa and Kcore variabilities on the spatial filtering. 6.1 Fi lter function for different specific storages In Chapter 2 we discussed that specific storage determines the influence zone, or, the filtering area of the slug test. In Chapter 5 we manifested the relationships Re ~ Ss and G(rD) ~ Ss by a perturbation method. In chapter 4 we presented a Wiener filtering methodology to determine the filter functions in heterogeneous fields. Now we utilize the Wiener filtering method and conduct a systematic investigation of the role of Ss on the filtering in heterogeneous fields. 134 Chapter 6. FILTER FUNCTION 135 6.1.1 Slug tests in isotropic core-scale fields: different field sizes, different specific storages We perform slug tests in an isotropic exponential K™"* field with Ss = 10 - n , n = 1, 2, . . . , 5. The field is generated by a sequential Gaussian simulation method of GSLIB [37] with a variogram -y(h) = c-[ 1 - e<~ »>], (6.120) where c is variance (c = 0y - c ) , and a is the correlation length in X and Y direction(a = Xx = Xy). The parameters are chosen to be, a\c = 0.25, Xx = Ay = 4Ax w , field size = 512 x 512 blocks(Aa:™ x Ax™).1 The cases of other parameter values or spatial variabilities are to be discussed in the next section. With different specific storages, slug tests are performed in this field at areas of 64 x 64, 128 x 128 and 256 x 256 blocks, respectively, listed in table 6.5. Sampling distance equal to 1 means slug test performed at every block in the selected area, equal to 2 means at every other blocks. Table 6.5: Slug tests with different specific storages. specific storage slug test area sampling distance 10"1, 0.05, 10"2, 0.005, 10"3, 10"4, 10"5 64 x 64 I 10"1, 10~2, 10"3, 10~4, IO"5 128 x 128 2 1Q-1, IO"2, 10"3 256 x 256 2 Fig. 6.48 shows the grey-scales and histograms of the core-scale and the measured field of size 64 x 64. The grey-scales clearly demonstrate the smoothing effect of the filtering process while this effect is embodied in the histograms as the shrink of data ranges. The smoothing is stronger for smaller specific storages due to a larger area of aquifer being averaged. 1 For illustration convenience, the units of correlation length and field blocks may be omitted in the rest of this Chapter. Chapter 6. FILTER FUNCTION 136 Figures 6.49, 6.50, 6.51 are the plots of the variograms of the core-scale and the measured fields of size 64 x 64, 128 x 128 and 256 x 256. We see that the averaging reduces the variance and increases the correlation range. The changes in mean and variance versus specific storage are plotted in Fig. 6.52. We see that the spatial filtering with different specific storages is approximately mean-preserving, i.e., the means of the measured fields(/iym) are close to the mean of the core-scale field(//yc). The deviation of pym from PYC may result from the the differences in the spatial filtering of each slug test excited in differing aquifer area of dissimilar local heterogeneity.. The variance decreases with decreasing Sa, because smaller Sa invokes a larger averaging area, i.e., a stronger 'smoothing effect'. Note that the means and variances of the measured fields with Sa = 10 - 4, 10 - 5 appear very similar. This is because the influence radii of slug tests with Sa < 10 - 4 have reached the whole inner heterogeneous zone(65 x 65 blocks) simulated by the slug test FDM model, thus their averaging effects become similar. Chapter 6. FILTER FUNCTION 137 Core-acale Log10(K) Manured Log10(K) Ss=lE-1 Measured Log10(K) Ss=1E-2 Measured Log10(K) 3»=1E-3 Measured Log10(K) Sa=1E-4 Measured Log10(K) Ss=1E-S Figure 6.48: The grey-scales and histograms of core-scale and measured fields with differ-ent specific storages. The core-scale field has an isotropic exponential variogram. Field size 64 x 64 blocks. Chapter 6. FILTER FUNCTION X-direction variogram of core-scale and measured fields 0.20 0.15 0.10 0.05 h 0.00 G B B B • J 18 20 Y-direction variogram of core-scale and measured fields 0.20 r 0.15 h 0.10 0.05 h 0.00 Figure 6.49: The variograms of core-scale and measured fields of size 64 x 64. Chapter 6. FILTER FUNCTION 139 Figure 6.50: The variograms of core-scale and measured fields. Field size 128 x 128. Chapter 6. FILTER FUNCTION 140 X-direction variogram of core-scale and measured fields Y-direction variogram of core-scale and measured fields 0 2 4 6 8 10 12 14 16 18 20 Distance Figure 6.51: The variograms of core-scale and measured fields. Field size 256 x 256. Chapter 6. FILTER FUNCTION 141 Figure 6.52: The mean and variance of core-scale and measured fields with different specific storages. Field size 64 x 64, 128 x 128, 256 x 256. Chapter 6. FILTER FUNCTION 142 6.1.2 Comparisons of Kcore and Kslu9 spectra In this section we analyze the spectra of Kcore and K*lug fields of different field size and different specific storages. The strength of spectral analysed is the visualization of the information contained in the data record in the frequency domain. The power spectrum density(PSD) illustrates the distribution of variance over spatial frequency [Gelhar, 1994, where R(0) is the autocovariance of zero lag. It also means that spectrum volume equals to the total energy of the data record -the summation of the square of the data at every point, divided by data record length. For a two dimensional periodogram, the PSD at zero frequency equals to the square of the volume of the 2D data series in physical space. In the following spectrum estimation by the multitaper method, the data mean is removed. The above properties may not be exactly met due to the way of data tapering and the removing of the mean. The estimated spectra may only demonstrate these properties approximately. The spectra of the core-scale fields and the measured fields with different specific storages are shown in Fig. 6.53. The graphs in Fig. 6.53 demonstrate: 1) The sharpening effects of the spectra for decreasing specific storages. In physical space it means the increasing of correlation range, because of the stronger averaging effect of smaller Ss. The first frequency shown is the Rayleigh frequency.2 2) The Spectra of S3 = IO - 4 , IO - 5 are very similar, this is due to the same reason for the similarity of the mean or variance between these measured fields, which we observed Fig. 6.52. The comparisons of the spectra of the core-scale and the measured fields of size 64 x 2 T h e zero frequency is removed from these graphs due to the logarithmic plotting requirement. This wil l be the case for the rest log scale P S D plots in this thesis. P31] *(0) (6.121) —oo Chapter 6. FILTER FUNCTION 143 64, 128 x 128, 256 x 256 blocks are shown in Fig. 6.54. The target spectrum is eq. 4.80 that is S(fx,fy) = 2 * X x V r . (6-122) Eq. 6.122 is the Fourier transform of the covariance function associated with the vari-ogram eq. 6.120, Cov(h) = c e<- (6.123) where c = cr2 and a = Xx = Xy. Thus the parameters of the target spectrum should be: a2 = 0.25, \ x = \ y = 4Axw. We see that the spectra of the core-scale fields deviate from the target spectrum at / > 0.2. This is caused by the inherent routines of GSLIB [37] for random field generation, not by the multitaper spectrum estimation method. We shall illustrate our explanation shortly. In Fig. 6.53 and Fig. 6.54 we also see that the spectra of the measured fields of size 64 x 64 have tails at / > 0.38 where energy drops off much less rapidly. Three possible causes for the tails are: 1) they are somehow transferred from the tail of the spectra of the core-scale fields; 2) there are some numerical errors in the measured fields, resulting from the slug test simulation model(Chapter 2) and the Yc inversion method (Chapter 3); and 3) they come with the multitaper method. These aspects are investigated through a series of hypothetical filtering tests. Firstly, filter the core-scale field of size 256 x 256 with the Gaussian filter of eq. 4.81 of filter widths Xx = Xy = 4, 8, 16, 32Axw, respectively. Secondly, select the filtered fields of size 64 x 64, 128 x 128 and 256 x 256,3 and estimate their spectra by the multitaper method. 3The fields of 64 x 64, 128 x 128 are chosen from the center area of the field of 256 x 256. Chapter 6. FILTER FUNCTION 144 Fig. 6.55 shows the spectra of core-scale fields and the hypothetically filtered fields by different filter widths. It displays that the spectra of the fields filtered with A = 4 have no tail; while the rest have tails. Consequently, 1) the foregoing observed tails are not necessarily affected by the tail of the spectra of the core-scale fields; 2) these tails do not have to be caused by numerical errors in the measured fields(i.e., these hypothetical filtering can be assumed containing negligible errors, but tails still observed.); 3) they are not likely coming with the spectrum estimation method(i.e., otherwise, one might observe tails consistently.). " ~ In Fig. 6.55, we also see that the spectral energy at the Rayleigh frequency, of different field sizes, has a very similar pattern as that in Fig. 6.53. This indicates that the nature of conductivity filtering of the slug test may be similar to a Gaussian filtering. From Fig. 6.55, we see that with increasing field size, the spectra of the filtered fields match better to the target spectra. This is evidently demonstrated in Fig. 6.56. From Fig. 6.56, we see that the pattern of spectral energy distribution at low frequency / < 0.07 is very similar to that of Fig. 6.53, even though spectral estimates have low reliability at low frequencies. Furthermore, a random field is generated with the spectral method introduced in Section 4.2.3.B on page 95. The parameters are chosen to be the same as the above core-scale field. Similarly, we perform the hypothetical filtering and estimate the spectra. Results are shown in Fig. 6.57, 6.58. The formats of Fig. 6.57, 6.58 are identical to those of Fig. 6.55, 6.56. Here we see that the PSD of the spectrally generated core-scale field matches the target spectrum very well. Other spectra characteristics are very similar to those of Fig. 6.55, 6.56. This indicates that the tails shown on the spectrum of the core-scale field generated by the GSLIB method should not be a problem and the core-scale field has no non-workable statistics. Chapter 6. FILTER FUNCTION 145 Now, the reason accounted for the spectral tails has to be related to the scale of filtering. Fields filtered with large filter width, equivalent to slug testing with a small Ss, will have corresponding large correlation lengths in physical space and sharp spectra in frequency space; consequently, it requires a fine resolution, i.e., a small Rayleigh frequency. Fig. 6.59 shows the Rayleigh frequencies and the target spectra of the above core-scale and filtered fields. On the other hand, these spectral tails have very small magnitude, in which case the log-log spectrum graphs at high frequency may be somehow misleading, for instance, the target spectra goes straightly down to ^ , whereas the spectrum volume(total energy) is finite, i.e., the energy above certain frequency should be no more than zero. Thus if the spectra is plotted in non-log space, one should observe the spectral energy approaching zero beyond certain frequency. Chapter 6. FILTER FUNCTION 146 spectrum of measured field size 256x256 IMOO h 1»04H 0.01 0.1 FREQUENCY spectrum of measured field size 128 x 128 ss = io- 5~io- 1 spectrum of measured field size 64 x64 S s = 10- 5 ~10 1 I M O O le-044 core-scale Ym:S»=le-l Ym:Ss=le-2 Ym:S»=lo-3 Ym:Ss=le-4 Ym:S»=le-5 0.01 0.1 FREQUENCY lfttOl 1C400 g leOl le-02 le-03 1 & 0 4 0.1 FREQUENCY Figure 6.53: Spectra of the core-scale fields and the measured fields with different specific storages. The core-scale field has an isotropic exponential variogram with al2%K = 0.25, Xx = \ y = AAxw. Chapter 6. FILTER FUNCTION 147 le+01 b-le4«0 t-spectra of core-scale fields lo+Ol h I04OO t-leOl le-02 1»04 spectra of measured fields ss = io-1 FREQUENCY FREQUENCY le+01 h le+OO g 1»01 le«2 le03 le04 spectra of measured fields Sfi = 10" 2 g 1*01 L lo03 h 1&04 spectra of measured fields Ss = 10- 3 FREQUENCY 0.01 0.1 FREQUENCY Figure 6.54: Comparisons of the spectra of core-scale fields and measured fields of size 64 x 64, 128 x 128, 256 x 256 blocks. Chapter 6. FILTER FUNCTION 148 spectra of core-scale and filtered fields core-scale: generated with GSLIB filtered: X - 8,16,32 core-scale filtered: X = 8 le+01 l f r O l t-le-02 t-le+00 1 * 0 2 J l e * H 1*0<H le-08-l le-10 h FREQUENCY FREQUENCY filtered: X = 16 filtered: X = 32 l e O U l * 0 4 j le06J l e O W le-10 le+OO k l e O Z H le-04-| 1*064 1*0*4 le-10 FREQUENCY 0.01 0.1 FREQUENCY Figure 6.55: Spectra of the core-scale fields and the hypothetically filtered fields with different filter widths. The core-scale fields are generated with a sequencial Gaussian method of GSLIB, with an exponential variogram, <rfgK = 0.25, Xx = Xy = 4Axw. Chapter 6. FILTER FUNCTION 149 field size: 256 x 256 spectra of core-scale and filtered fields core-scale: generated with GSLIB filtered: X = 4,8,16,32 1*400 le<ttJ le-04- | le-06-l l e - 1 0 F R E Q U E N C Y field size: 128x128 field size: 64 x 64 le -024 le-04-| le -064 l e - 1 0 r-• TargetO Target4 Targets Ta rge t l6 . Target 32 Core - l ca l e - • — 1 j t m t a = 4 -i 1 a m r f w = g -w— L a m d a = 16 - * \ \ a 1 . \ L a m d a = 32 1\ 1 1 \ i u \ \ \ Xx v a. \ V Vv l * > » * L ft \ \ k \ \ \ ^^NSJL * 11 \ V 0.01 0.1 F R E Q U E N C Y le+OO I-1&02J 1&04-I l e - 0 6 j l e - 0 8 J l e - 1 0 h TargetO ^ \ \ / V » »<. T a r g « 4 \ \ \ \ Targets \ \Y», Targe t l6 \ TVS T a r g « 3 2 \ \ \ Core^icale — \ \ Lamda = 4 - t— \ \ L a m d a = 8 — \ \ \ \ ** V L a m d a = 16 - * — > • L a m d a = 32 - • — \ \ \ \ \ \ \ \ 1 \ 1 •• 1 1 -0.01 F R E Q U E N C Y 0.1 Figure 6.56: Comparisons of the spectra of the core-scale field generated with GSLIB and the hypothetically filtered fields of size 64 x 64, 128 x 128, 256 x 256 blocks. Chapter 6. FILTER FUNCTION 150 spectra of core-scale and filtered fields core-scale: generated with a spectral method filtered: X = 8,16,32 core-scale filtered: X = 8 0.01 0.1 0.01 0.1 FREQUENCY FREQUENCY filtered: X=16 filtered: X = 32 0.01 0.1 0.01 0.1 FREQUENCY FREQUENCY Figure 6.57: Spectra of the core-scale fields and the hypothetically filtered fields with different filter widths. The core-scale field is generated with a spectral method with an exponential variogram, cr?K = 0.25, \ x = \ y = 4Axw. Chapter 6. FILTER FUNCTION 151 field size: 256 x 256 spectra of core-scale and filtered fields core-scale: generated with a spectral method filtered: X = 4,8,16,32 le+OO h le02j 1&04-I 1&064 le-08-l l e - 1 0 TargetO Target* Torget8 T a r g e U 6 Target32 Core-fcale I mnAn — 4 8 L a m d a " 16 L a m d a = 32 FREQUENCY field size: 128 x 128 field size: 64 x 64 le+OO h 1 & 0 2 J 2 le*U le-06-l le-08-l l e - 1 0 TargetO T a r g e U Targets Target 16 Target32 Core-scale * — = 4 -+— L a m d a = 8 -**— L a m d a = 16 - * — L a m d a = 32 le+OO l e - 0 2 J 5 le-044 1&06H le -08 l e - 1 0 TargetO Target4 Targets T a i g e t l 6 Target32 Core-scale L a m d a = 4 I jimHw s 8 L a m d a = 16 L a m d a = 32 0.01 FREQUENCY 0.1 0.01 FREQUENCY 0.1 Figure 6.58: Comparisons of the spectra of the spectrally generated core-scale and the hypothetically filtered fields of size 64 x 64, 128 x 128, 256 x 256 blocks. Chapter 6. FILTER FUNCTION 152 target spectra of core-scale and filtered fields l e + 0 4 r . . — • . — . . — le-03 le-02 le-01 le+OO FREQUENCY Figure 6.59: The Rayleigh freuencies and the target spectra of core-scale and filtered fields. Chapter 6. FILTER FUNCTION 153 6.1.3 Comparisons of filter spectra, filter functions Fig. 6.60 shows the filter spectra with different specific storages. This time the frequency-coordinate is in non-log form. We can observe the sharpening effects of the spectra for decreasing specific storages, in physical space meaning that the filtering occurs in a large influence zone. Note that the spectral energy at zero frequency equals to the square of the filter volume because the filter spectrum is estimated as a periodogram in the iterative Winer filtering approach. We see that this energy decreases with Sa. This is because smaller Ss introduces lower filtering power and a larger influence zone, however, the effective filtering is limited to the 65 x 65 blocks in the inner heterogeneous zone of the simulation model. Fig. 6.61 shows the comparisons of the filter spectra of the fields of size 64 x 64, 128 x 128, 256 x 256. Their differences at low frequency were caused by the role of the Rayleigh frequency in the spectrum estimation, which is different for the measured fields that have different correlation lengths. We see that the filter spectra of the measured fields have a good match for Ss = 0.1. This means the Rayleigh frequency ^ is good enough for the spectrum resolution of the measured fields with Ss = 0.1. For Ss = 0.001, we see that even the Rayleigh frequency still appears to be too coarse, which is evidenced by the fact that the energy at zero frequency is less than 1. The filter function comparisons are illustrated in Fig. 6.62. The first three graphs show the comparisons of the filter functions associated with the slug tests of size 64 x 64, 128 x 128 and 256 x 256 with the filter functions identified by the perturbation method in Chapter 4. The filter functions by the spectral method are radially averaged. For a fair comparison to the spectral results that only have discrete values at dimensionless radius(scaled by wellbore radius) 4 , rD = 2, 4, • • • 64, the filter functions from the 4 F r o m now on, for the convenience for plotting and description, the wellbore radius is approximated by half of a block size, it actually, in the slug test simulation model, equals to 0.63Ax, 0.599Aa;, 0.582Ax Chapter 6. FILTER FUNCTION 154 perturbation method are averaged as / G(s)ds G(rD) = ^ 2 — (6.124) Note that the filter functions from these two methods are in fact not accurately compara-ble. The perturbation results are obtained from the small perturbations in a homogeneous aquifer, thus the flow pattern in the perturbed aquifer is very close to that of a homo-geneous case, whereas the filter identified by the spectral method can be interpreted as the one averaged from all individual filter functions of every slug test performed in the heterogeneous core-scale field. Also note that the radius rD in the filter function of the spectral method are greater than its true value(see the footnote of page 153), the filter functions would match more closely to the perturbation results if the accurate rD is used. The first three graphs show that the filter functions of the fields of size 128 x 128 and 256 x 256 match very well and they are closer to the filter functions of the perturbation method than are the filter functions of the field of size 64 x 64. The last three graphs display the filter function comparisons for different storages. Note that their differences can not be well shown on the scale of a same plot. For instance, the influence radii Re or the equivalent filter widths Ae appear to be very similar, whereas from the perturbation results we see that they are very different, e.g., Re = 19.5, 60, 195, Ae = 5, 13, 20, for sa = i o - \ i o - 2 , i o - 3 . Fig. 6.63 shows the filter functions in space where the central point corresponds to the location of the wellbore. The filtering power is magnified by 104. The 20-contour demarks an area where the filter is well resolved. We see that filter function of a smaller Sa has less spurious off-center peaks, because filter is better resolved for a smoother field. Also, the the filter in the field of size 256 x 256 has less off center peaks than that of the field for S,=0.1, 0.01, 0.001, respectively. See Fig . 2.9 Chapter 6. FILTER FUNCTION 155 of size 128 x 128. Because the filter identified by the spectral method can be interpreted as the one averaged from all individual filter functions of every slug test performed in the core-scale field, thus, the more slug tests are performed, the more averaging the spectral method does to the filter function. The corresponding filter functions by the perturbation method are shown in Fig. 6.64. It can be seen from Fig. 6.64 and Fig. 6.63 that the 20 contour line of S3 = 10 - 2 is larger than that of Ss = 10 _ 1 while the 20 contour of Ss = 10~2 and S, = 10 - 3 stay approximately the same. For Ss decreasing from 10 - 1 to 10 - 3, the 220 contour line shrinks. Filter functions from the perturbation method clearly demonstrate the differences of the averaging power less than 1 for different specific storages while such differences become very spurious in the spectral method. This existence of such an averaged filter identified by the spectral method can be directly evaluated through the scatter cross-plots between the measured field, Ym, and the estimated field from Winer filtering, Ym, which is the result of the core-scale field filtered by the optimally estimated filter function(see Fig. 6.65, 6.66, 6.67). If these scatter points gather around the line of Ym = Ym, a universal filter for all the individual slug test exists. Note that the measured Ym contains some noises from the simulation model and the inversion procedure. These noises do not form up to a trend, for instance, in the inversion procedure different Ym estimates can fit the requirement of the objective function. Thus the cross-plots of Ym versus Ym need not to be exactly on the Ym = Ym line. If the cross-plots obviously skewed away from the Ym = Ym line, a universal filter does not exist (the filter is poorly resolved) or the filtering is nonlinear, or, there are some trends in the scale-scale and measured Yc data. This may be the case where the core-scale field has a large variance or exhibits a multimodal distribution, e.g., the bimodal sand-shale lense field. Fig. 6.68 displays the comparison of the equivalent filter width versus specific storage Chapter 6. FILTER FUNCTION 156 (Ae ~ is(^sl)) fr°m the spectral analysis and the perturbation method. In the perturba-tion method the equivalent filter width is defined in eq. (5.110) which is A« = wfir (-^r <6-125> In the spectral analysis method, the filter volume Vo(Re) is directly obtained by the value of transfer function at zero frequency, i.e., Va(Re) = G(0). Va(Re) is found to be close to 1, meaning that the filtering process is approximately mean-preserving. G(l + £) equals G(2), because the closest node to the center of the wellbore is located at dimensionless distance rD = 2rw/rw = 2. This explains why the equivalent filter widths of the spectral method are greater than that of the perturbation method.5 We also see that the Ae of Spectral(128 x 128)6 is greater than that of Spectral(64 x 64), because, in calculating G(0), Spectral(128 x 128) has a finer Rayleigh frequency, which promotes the filter volume, i.e., G(0). Furthermore, we see that the Ae of Spectral(128 x 128) and Spectral(256 x 256) become similar, this reflects the fact that refining the Rayleigh frequency becomes insignificant in promoting G(0). On the other hand, this tells us that the Ae of Spectral(256 x 256) can be considered as the one of the highest accu-racy. We fit the Ae of Spectral(256 x 256) versus S„, i.e., Ae ~ lg(^j), by a parametric representation, Ae = A° - [ l g (^=) ] ° + 14-58 l g ( J L ) (6.126) A" = 12.5, [lg(-^=)]° = 0.5. Note that this relationship is obtained through the slug tests in the fields with variance <rYc = 0.25, and correlation length A x = Aj, = 4AZ*". 5Because in the perturbation method in Chapter 5 we have observed that VG(RC) « 1 and G(l + £) « G(l), f o r £ ^ 0 + . 6'SpectraP denotes thatAe is identified by the spectral analysis method, '128 x 128' represents the slug testing field size. Chapter 6. FILTER FUNCTION 157 In the following, we present a systematic exploration of the spatial variability of Yc on the spatial filtering. Chapter 6. FILTER FUNCTION 158 filter spectrum of loguK field of size 128 x 128 ss = io-5~io-1 T 1 1 1 r I I : I I I i _ 0 0.1 02 0.3 0 4 0 5 F R E Q U E N C Y filter spectrum of log^K field of size 64 x 64 ss = io-5-io-1 i 1 i r -I 1 I I I . L 0 0.1 02 0.3 0 4 0 : F R E Q U E N C Y Figure 6.60: Filter spectra of the slug tests in the conductivity fields of different specific storages. Chapter 6. FILTER FUNCTION 159 comparison of filter spectra of slug tests in fields of different size 0.8 h 0.6 h OA k filter spectrum 88 = 10-! filter spectrum s s = io- 1 0.8 h 0.6 r-0.4 h 0.2 h filter spectrum s s = io- 3 0.1 02 0.3 F R E Q U E N C Y Figure 6.61: Comparisons of the filter spectra in fields of size 64 x 64, 128 x 128, 256 x 256 blocks. Chapter 6. FILTER FUNCTION continued on next page • Chapter 6. FILTER FUNCTION Figure 6.62: Averaged filter functions for different specific storges. Chapter 6. FILTER FUNCTION 162 filter function Ss = 1Q-* [256x256] filter function Ss = 10-» [128x128] 0 c 1 < O < o A -> \ 1 / 0 A 1 1 V W) U V 0 > < J < o O n V V < 7^ o cz > g o < c : V O /I 1 j 0 V /))}\ o V 0 C ? < D < o o 0 o „..,• -,-tO> 0 -aa filter function Ss = lQ-2 [256x256] > n <• 0 r • (Q A rs 0 *. \ 4£ 0 0 s c ) filter function Ss = 10-2 [128x128] O < • V - <=> < o O > 0 0 Ct 0 0 0 V! J o c 0 o • 0 0 o filter function Ss = IO 3 [256x256] 0 r I J s Z^ ). filter function Ss = 10-3 [128x128] Figure 6.63: Filter functions for different specific storges from spectral method. Enlarged by 104. Contour level = 200. Chapter 6. FILTER FUNCTION 163 Ss = io-i Figure 6.64: Filter functions for different specific storges from perturbation method. Enlarged by 1G4. Contour level of inner zone = 200; contour level of outer zone = 0.2. Chapter 6. FILTER FUNCTION 164 Figure 6.65: Scatter plots of measured fields versus estimated fields from Wiener filtering. Field size 256 x 256. The inclined lines are of unit slope. Chapter 6. FILTER FUNCTION 165 Figure 6.66: Scatter plots of measured fields versus estimated fields from Wiener filtering. Field size 128 x 128. The inclined lines are of unit slope. Chapter 6. FILTER FUNCTION 166 Figure 6.67: Scatter plots of measured fields versus estimated fields from Wiener filtering. Field size 64 x 64. The inclined lines are of unit slope. Chapter 6. FILTER FUNCTION 167 80 0 i— •—. — , — , • 1 — le+01 le+02 le+03 le+04 le+05 Specific storage[l/sqrt(Ss)] Figure 6.68: Comparison of the equivalent filter width versus specific storage (Ae ~ lg(^j)) from the spectral analysis and the perturbation method. Chapter 6. FILTER FUNCTION 168 6.2 Fi l ter functions in conductivity fields of different spatial variabilities The statistical parameters describing the spatial variability of a conductivity field may include the mean, variance, correlation length, variogram model, and stochastic realiza-tion. By the scaling property K® and the spatial filtering representation of eq. (3.29), the mean should not affect the filtering, thus it needs not be considered. The other effects of the spatial variability on the conductivity filtering are discussed in the following sections. 6.2.1 Variance, variogram model and stochastic realization We investigate the role of the variance, variogram and stochastic realization in the filtering process by the slug tests in isotropic conductivity fields(l^) of multiple realizations. 6.2.1. A Mult ip le realizations of isotropic fields The multiple realizations of isotropic fields are generated by a sequential Gaussian simu-lation method of GSLIB [37] with both the exponential variogram (6.120) and a Gaussian variogram [37], >y(h) = c- [ 1 - e ( _ S>], (6.127) The parameters are chosen to be a = \ x = \ y = 8Axw, c = o~Yc — 0-3, 0.7, field size = 256 x 256 blocks. With Sa = 0.001, sampling distance = 2, we perform slug tests in a center area of 128 x 128 blocks in these fields.(Table 6.6). Fig. 6.69 shows the estimated filter functions in the above fields. The last graph is the plot of the respective averaged filter functions from the first three graphs. Their close-match indicates that the variance, variogram model and multiple realization only have very small effects on the filtering process. This shows their influence on the changes of the filtering pattern of individual slug tests are averaged out by the spectral method. Chapters. FILTER FUNCTION 169 Table 6.6: Slug tests in isotropic core-scale fields from multiple realizations. variogram model variance(<Ty) number of realization exponential 0.3 3 exponential 0.7 3 Gaussian 0.3 3 Thus the universal filter from spectral analysis is approximately independent of variance, variogram model and statistical realization. Note that this observation is obtained from isotropic fields with modest variance a\c < 0.7. The role of anisotropy will be discussed in the next section. Chapter 6. FILTER FUNCTION filter function exponential fields <*2LoglO(K) = 0-3 filter function exponential fields ° 2 LoglO(K) = u - 7 filter function Gaussian fields ° 2 LoglO(K) = O-3 E x p o n e n t i a l f i e l d l [ V a r i a n c e = 0 . 3 ] - • — E x p o n e n t i a l field 2 [ V a r i a n c e = 0 . 3 ] - t — E x p o n e n t i a l field 3 [ V a r i a n c e = 0 . 3 ] - a - -20 E x p o n e n t i a l field l [ V a r i a n c e = 0 . 7 ] -e— E x p o n e n t i a l f i e l d 2 [Var i ance=0 .7 ] -+— E x p o n e n t i a l field 3 [Var i ance=0 .7 ] - a - -20 G a u s s i a n field l [ V a r i a n c e = 0 . 3 ] - • — G a u s s i a n field l [ V a r i a n c e = 0 . 3 ] G a u s s i a n field l [ V a r i a n c e = 0 . 3 ] - a - -D i m e n s i o n l e s s d is tance 2 0 continued on next page • Chapter 6. FILTER FUNCTION 171 0.04 averaged filter functions of the exponential and Gaussian fields 5 10 15 20 Dimensionless distance Figure 6.69: Filter functions estimated in multiple core-scale conductivity realizations, with different variances and variogram models. The last graph is a comparison of the averaged filter functions shown in the first three graphs. Chapter 6. FILTER FUNCTION 172 6.2.2 Core-scale conductivity of different correlation length We explore the role of the correlation length in the filtering with two slug test schemes. One is slug tests in isotropic fields with different correlation lengths, the second is slug tests in anisotropic fields with different correlation length ratios. 6.2.2.A Isotropic, different correlation length The isotropic fields are generated by a sequential Gaussian simulation method of GSLIB [37] with the exponential variogram of equation (6.120). The parameters are chosen to be a = \ x = Xy = 3, 6, 12Axw, c = <rYc = 0.3, field size = 256 x 256 blocks. With Sa = 0.001, sampling distance = 2, we perform slug tests in a center area of 128 x 128 blocks. Fig. 6.70 shows the comparisons of the estimated filter functions. Their close-match demonstrates that the power distribution of the Yc filtering is approximately independent of the correlation length of isotropic fields. Fig. 6.71 shows the scatter plots of measured fields versus estimated fields from Wiener filtering. Here we see that the filter is better resolved with increasing correlation length. This reflects the fact that the similarity of the spatial filtering among individual slug tests increases with correlation length. Chapter 6. FILTER FUNCTION 173 comparions of filter functions in isotropic IgK fields of different correlation lengths X 0.06 0.05 I-0.04 0.03 0.02 0.01 0.00 L A M B D A = 3 - H — L A M B D A = 6 -o— L A M B D A B 1 2 10 IS DIMENSIONLESS DISTANCE 20 Figure 6.70: Comparions of filter functions in isotropic fields of different correlation lengths, Xx = Xy = 3, 6, 12Axw. Chapter 6. FILTER FUNCTION 174 A = 12 scatter plots of measured fields versus estimated fields from Wiener filtering core-scale IgK fields: isotropic, different correlation X Ss = 0.01 - 2 - 1 0 1 2 Measured Ym X = 6 X = 3 Measured Ym Measured Ym Figure 6.71: Scatter plots of the measured fields versus the estimated fields from Wiener filtering. The core-scale fields have isotropic vaiograms of correlation lengths A = 3, 6, 12Axw. Chapter 6. FILTER FUNCTION 175 6.2.2.B Anisotropic, different correlation length ratios The anisotropic fields are generated by a sequential Gaussian simulation method of GSLIB [37] with the exponential variogram (6.120). The parameters are chosen to be: <rYc = 0.3, Ay = 3Axw, field size = 256 x 256 blocks. The correlation length ratios are listed in Table 6.7. With Sa = 0.01 and sampling distance = 1, slug tests are performed in a center area of 64 x 64 blocks in these fields. Table 6.7: Slug tests in anisotropic conductivity fields. fields of different Xx/Xy ratios \ y = ZAxw Xx/Xy 1 3 5 10 20 30 40 Fig. 6.72 shows the grey-scales and histograms of the core-scale and the measured fields. We see that the averaging effect in the X direction is stronger than that in the Y direction, and it is more pronounced for increasing Xx/Xy ratio. Fig. 6.73 shows the reduced variance, aYc — o~Ym-> v e r s u s the correlation length ratio It demonstrates that the stronger the anisotropy, the less the variance reduction. Because the anisotropy somehow bars the flux to the wellbore in the Y direction, and forces the averaging to occur in a relatively long and narrow zone, so that the net averaging volume may be relatively smaller and variance is less reduced. Fig. 6.74 shows filter functions on the the center-lines in the X and Y directions. We see that X-direction averaging expands to a larger area and y-direction averaging stays approximately the same. This phenomenon is shown in space in Fig. 6.75. The scatter plots in Fig. 6.76 indicates the existence of a universal filter for the individual filtering of the slug tests performed in the corresponding fields. Note that all the Yc fields are generated with a fixed Ay = 3Ax*", if we increase A y but keep the Chapter 6. FILTER FUNCTION 176 same Xx/Xy ratios, such averaged filters identified with the spectral method may not be meaningful, i.e., the averaging patterns of slug tests performed at different locations may be very dissimilar, hence, degrades the spectral results. Chapter 6. FILTER FUNCTION continued on next page • Chapter 6. FILTER FUNCTION 178 Figure 6.72: The grey-scales and histograms of the core-scale and the measured field of different correlation length ratios. The fields shown are 64 x 64 blocks each. Chapter 6. FILTER FUNCTION 179 1 9 4> s I eg & I I V u a .a es correlation length ratio Xx/Xy 1 3 a •s a es a -0.02 -0.04 -0.06 -0.08 -0.10 I--0.12 -0.1 -0.08 -0.06 mean of core-scale -0.04 Figure 6.73: Graph A: reduced variance(cJyc — oYm) versus correlation lenght ratio Graph B: the mean of measured fields pYm versus the mean of core-scale fields pYc. Chapter 6. FILTER FUNCTION 180 0.16 0.14 0.12 0.10 | 0.08 i filter function in isotropic field X, = Xy = 6 r w X center-line Y center-line l l -16 0 16 Dimensionless distance 32 I filter function in ansotropic field X,/Xy = 3,Xy = 6r„ X center-line Y center-line • -16 0 Dimensionless c 0.16 0.14 0.12 g " i o 0.08 0.06 0.04 0.02 0.00 text* filter function in ansotropic field Xx/Xy = 5, Xy = 6 r w X center-line ' Y center-line ' -32 filter function in ansotropic field Xji/Xy = 30, Xy = 6r„ -16 0 Dimensionless • t i« t » n<y -16 0 Dinuns ion less distance Figure 6.74: The filter functions on the center lines in the X and Y directions in anisotropic conductivity fields. Chapter 6. FILTER FUNCTION 181 filter function in ansotropic field X x / X y = 5,Xy = 6 r w -o o • 0 v < o \ •—t— 0> o O V • 0 0 ^ < > <\) 6 o o 32 -32 32 32 i filter function in ansotropic field A,/Xy = 30,Xy = 6r„ > < => < > > o o V o C I > • c— • * — -—-^ ^ c <= 3 Ofr Figure 6.75: Filter functions of different anisotropy ratios. Enlarged by 104. Chapter 6. FILTER FUNCTION 182 scatter plots of measured fields versus estimated fields from Wiener filtering isotropic field X» = Xy = 6 r w anisotropic field X,/Xy = 3,Xy = 6r, -05 0 05 Measured Ym -05 0 05 Measured Ym anisotropic field Xj / Xy = 5, Xy = 6 r w anisotropic field Xx/Xy = 30,Xy = 6r„ -05 0 05 Measured Ym -05 0 Measured Ym Figure 6.76: Scatter plots of the measured fields versus the estimated fields from Wiener filtering. Anisotropy ratio of the core-scale fields: Xx/Xy = 1, 3, 5, 30Axw. The inclined lines represent unit slope. Chapter 6. FILTER FUNCTION 183 6.2.2.C Lense field We can generate a bimodal lense field by several steps: Firstly, generating two isotropic fields, denoted as YCA, YB, with a sequential Gaussian simulation method of GSLIB [37]; secondly, generating the lense shapes with a Boolean simulation method [37]; lastly reducing the mean of YCA and projecting it through the Boolean lense shapes to YB. Two lense fields are generated by this approach, see Table 6.8. Slug tests are performed in the two lense fields with different specific storages, see Table 6.9. Table 6.8: Statistic parameters for two lense fields. variograms of correlation of reduced mean of variance of YCA, YCB YC\ YCH YCA YC\ YCB field one exponential 8 -2 0.3 field two Gaussian 4 -4 0.8 Table 6.9: Slug tests in two lense fields. S8 slug testing area sampling distance field one IO"1 128 x 128 2 10~2 128 x 128 2 field two io- 1 64 x 64 1 IO"2 64 x 64 1 Fig. 6.77 shows that the grey-scales and histograms of the core-scale and the mea-sured fields. We see that the filtering occurs to both the background aquifer(sands) and foreground lenses(clay lenses), but there are only limited averaging between the sands and the lenses. Fig. 6.78 displays the estimated filter functions. We discussed in the previous sections that the filter function from the spectral method can be interpreted as an averaged filter. Chapter 6. FILTER FUNCTION 184 Clearly the filtering for slug tests in the lenses and those in the sands are fundamentally different. Thus, the averaged filter function by the spectral method is not representative to the individual filtering patterns of each slug test and it is poorly resolved. Fig. 6.79 demonstrates the scatter plots of the measured fields versus the estimated fields from Wiener filtering. It can be seen that the optimally estimated filter functions honor the filtering pattern in the sand aquifer. The pattern of the scatter plots reflects the principle filtering fashion, which is affected by the extent of the lense heterogeneity. The 128 x 128 field is relatively less heterogeneous than the 64 x 64 field. The 64 x 64 field basically comprises half lenses and half sands, resulting in two distinctive centers in the scatter plot. Chapter 6. FILTER FUNCTION 185 Core-scale ler.se one [128x128] Measured Log10(K) Ss = 0.1 Measured Log10(K) Sa . 0.01 imrmu U . 1 0 K I uaiaw Figure 6.77: The grey-scales and histograms of the core-scale and the measured field with lense heterogeneties. Size of Lense Field one = 128 x 128; size of Lense Field two = 64 x 64. Chapter 6. FILTER FUNCTION 186 filter function in lense field one size: 128 x 128 Ss = 0.1 0 | % Si = 0.01 32 32 v / o (—\ 0 1 £ ^ 0 ^—Y\ - 2 0 - ^ Y ^ A 0 > .0 r i I/-? 0 -filter function in lense field two size: 64 x 64 Ss = 0.05 > *—.. 1 1—/ 3 (? . ^ -ft* Ss = 0.005 32 -32 Figure 6.78: Poorly resolved filter functions in two lense fields. Enlarged by 104. Chapter 6. FILTER FUNCTION 187 scatter plots of measured fields versus estimated fields from Wiener filtering lense field one [128x128] Ss = 0.1 Ss = 0.01 Measured Ym Measured Ym Figure 6.79: Scatter plots of the measured fields versus the estimated fields from Wiener filtering. The inclined lines are of unit slope. Chapter 7 C O N C L U S I O N S 7.1 T h e spatial filtering from core-scale conductivity to the conductivity measured by a slug test 1. System and filtering. We conceptualize the scale relationship between the core-scale conductivity and the conductivity measured by a slug test in a system and filtering framework. The slug test is a stable, space-variant, and nonlinear filtering system. The extent of space-variant property and nonlinearity increases with increasing <Tyc(variance of core-scale field), and increases with decreasing .^(specific storage). 2. System outputs - inverse parameter estimation, (lg Kh — lg Ka) ~ (lg Sb — lg Slrue) are positively correlated due to the role of lg Sb in controlling the shape of the fitting curve to match the response curve measured in a slug test. Superscripts ° and 6 denote the inverse estimation methods which estimate K alone and concurrently estimate K and S„, respectively. 3. Role of Ss in filtering. In heterogeneous media1, the spatial filtering, from core-scale conductivity to the conductivity measured by a slug test, is dependent on the specific storage. With decreasing S3, the filtering area increases and the filter am-plitude proximal to wellbore decreases, and the filtering power distribution, from wellbore to influence radius, decays off less rapidly. Such filtering processes are ap-proximately mean-preserving, and it reduces variance and increases the correlation 1 Here "media" refers to core-scale conductivity fields. 188 Chapter 7. CONCLUSIONS 189 length of the measured conductivity fields. The changes in the variance and in the correlation length is proportional to decreasing S3. For the slug tests in the fields with variance oYc = 0.25, and correlation length A x = A y = 4Axtu(wellbore block diameter), it is found that equivalent filter width versus S8, Ae ~ lg(^;), Ae = A0 - [ l g ( ^ = ) ] ° + 14.58 lg( J L ) (7.128) A0 = 12.5, [lg(-L)]o = 0.5. 4. Effects of spatial variability on filtering. In isotropic heterogeneous media, in tests with o~Yc — 0-7> Gaussian and exponential variogram, correlation length 3Axw < A < 12AX*", and a number of realizations for each set of parameters, it is found that spatial filtering is approximately independent of variance, variogram model and stochastical realizations; it is also independent of the correlation length, but, for fields with larger correlation, the filter function is better resolved in the spectral analysis method. In anisotropic heterogeneous media, the spatial filtering shows an anisotropic pattern proportional to the anisotropy ratio. 5. Role of Ss in filtering. In homogeneous media, the spatial filtering is dependent on Ss- The filtering function has the form, G(rD) = G(1+0TAP, £ - 0 + , rD=r/rw (7.129) \rD) where the (?(!+£) is the filter amplitude at the wellbore edge, and p is an exponent: Ss io- 1 IO"2 io- 3 10 - 4 10"5 io- 6 IO"7 10"8 IO"9 1 0 - i o G(l+0 x IO"2 18.338 7.908 4.882 3.472 2.840 2.412 2.165 1.950 1.885 1.753 P 3.0 2.3 2.1 2.05 2.0 2.0 2.0 2.0 2.0 2.0 Chapter 7. CONCLUSIONS 190 where: 1) G(l + £) decreases with decreasing Sa, and 2) G(rD) ~ rD approaches approximately a l / r j decreasing law for Sa smaller than 10~4. The direct relation-ship between G(rD) and Sa is _ G(rD) = (0.2ff-25 + 0.2S?-75 + 0.01753) J ^ , , . (7-130) Equivalent filter width versus Sa, A e ~ lg(^-), Ae = A° - [ l g ( ^ = ) ] ° + 14.58 l g ( J L ) (7.131) A 0 = 5.47, [lg( J = ) ] ° = 0.5 Influence radius versus Sa, Re ~ l g(^r-), lgi?e = (lgi? e)° - [ l g ( ^ ) ] ° + (7-132) R e = 1 1 0 0»«-)° - W ^ f f H " ( 7 . 1 3 3 ) (lgi*er = 1.29, [lg(-^ )]<> = 0.5 7.2 Possible engineering application The spatial filtering methodology is equally applicable to many other engineering mea-surement techniques which invoke the spatial averaging of the measured properties. The measurement niters quantify the spatial significance of the sub-scale properties on the measured-scale properties. Here, we illustrate this idea in two types of problems in the following. 7.2.1 The Q-problem The Q-problem refers to the one whose primary concern is the flux rate under certain conditions. The scale under consideration in the Q-problem is the scale of the aquifer Chapter 7. CONCLUSIONS 191 that contributes the fluxes. The filtering property tells us that the measured conductivity (equivalent or effective conductivity) is strongly affected by near-wellbore heterogeneity, which thus is very crucial in the control of flux rate. Then direct implications may be in the siting of pumping wells or any other engineering facilities extracting water from aquifers. On the other hand, it is no wonder that the equivalent or effective conductivity estimated with traditional geostatistical methods, e.g., Kriging, conditional simulation, etc., which usually do not incorporate the filters invoked by engineering activities, might be far different from the one estimated with the filtering function, and thus result in errors in flux rate prediction. For instance, the statistically interpolated grid-scale conductivity may incur large error in the prediction of pump yield. 7.2.2 The X-problem The X-problem refers to the one whose primary concern is the mass transport distance under certain conditions. In the .//-problem the small-scale(or core-scale) conductivity is of practical interests, because the transport distance is often seen to be strongly im-pacted by the small-scale conductivity variations. The small-scale conductivity can not be measured or at least are too expensive to be measured, but if the filter function(from the small-scale to the measured scale) is known, one can deconvolve the measured data back to the small-scale using spectral analysis.2 On the other hand, when the measurements are insufficient, one may attempt to em-ploy some geostatistical models to simulate the the parameter values in the unsampled region, conditioning on the measured data. Since we now know from measurement filter-ing that the measured conductivity is a spatial average of the small-scale conductivity, 2 T h e result of the deconvolution is a spectrum or a variogram which can be used to generate small-scale fields. Chapter 7. CONCLUSIONS 192 thus the simulation does not have to: 1) exactly utilize the variogram inferred from mea-sured data record, or 2) exactly honor the measured data at the conditioning locations. In this case, deconvolution(to small-scale or to remove noise), or some soft-data techniques appear necessary. The soft-data techniques may incorporate all the available informa-tion of different characteristics and of various sources(like hard data and experience), thus may be more appropriate to relieve some engineering intractability. 7.3 Theoretical significance As we discussed in the introduction in Chapter 1, measurement filter is not only able to transfer parameter from small-scale to larger scale, but it can be employed in upscaling the fundamental flow and transport equations. The equations are often obtained from small-scale laboratory experiments and will become invalid beyond certain scale ranges. The measurement filter can be used to enlarge or transfer the application scale ranges of the equations, to accommodate the scale of field problems. 7.4 Future work 1. In the spatial filtering of heterogeneous core-scale conductivity, for each slug test we have only simulated an inner zone of 65 x 65 blocks(Axtw x Axw), and Ss = IO -", n = 1, 2, . . . , 5, due to the computational limitation. The simulation area and Ss value can be expanded if superior computer resources becomes available, or, more efficient numerical techniques(e.g., multigrid method) are developed. In unconfined aquifer where vertical flow becomes significant, 3D model should be used. 2. The filter function has not yet been parameterized in terms of the variance, corre-lation length and anisotropy of core-scale field. A generic form of filter functions possibly can be obtained if more tests could be performed with more computational Chapter 7. CONCLUSIONS 193 resource. So at this time, the nonparametric representation may not be readily ap-plied by field workers. But the generic formulas in homogeneous fields obtained with the perturbation method should serve as a good clue for the heterogeneous case. On the other hand, the similarity between the slug test filtering and Gaussian filtering, G,,,,) = « « - [ W + <*>" 1 (7.134) 7T Ax Ay where Xx and Ay are the characteristic filter widths which are somehow related to the variance, correlation length, anisotropy, and specific storage of the core-scale field, may ease the application dearly. 3. Only slug test filter is considered in the thesis. Other aquifer tests or engineering measurement techniques which invoke the spatial averaging of the intended-to-be-measured properties are worth being investigated and made readily applicable to field workers and theoreticians. Bibliography [1] Anderson, T. B., and R. Jackson, Fluid mechanical description of fluidized beds, Equation of motion, Ind. Eng. Chem. Fundam., 6(4), 527-539, 1967. [2] Ababou, R. and E. F. Wood, Comment on "Effective groundwater model parame-ter values: Influence of spatial variability of hydraulic conductivity, leakance, and recharge" by J. J. Gomez-Hernandez and S. M. Gorelick, Water Resour. Res., 26(8), 1843-1846, 1990. [3] Bachu, S. and D. Guthiell, Effects of core-scale heterogeneity on steady state and transient fluid in porous media: Numerical analysis. Water Resour. Res., 26(5), 863-874, 1990. [4] Bakr, A. A., Gelhar L. W., Gutjahr A. L., and J. R. MacMillan, Stochastical analysis of spatial variability in subsurface flows 1. Comparison of one- and three-dimensional flows. Water Resour. Res., 14(2), 263-271, 1978. [5] Barker, J. A., and R. Herbert, Pumping tests in patchy aquifers. Ground Water, 20(2), 150-155, 1982. [6] Baveye, P. and G. Sposito, The operational significance of the continuum hypothesis in the theory of water movement through soils and aquifers. Water Resour. Res., 20(5), 521-530, 1984. [7] Bear, J., Dynamics of Fluids in porous media, American Elsevier, New York, 1972. [8] Beckie, R. and B. Wang, A numerical method to characterize the averaging process invoked by a slug test, Computational Methods in Water Resources X, A. Peters et al. (eds.), 703-710, 1994. [9] Beven, K., Estimating transport parameters at the grid scale: On the value of a single measurement. J. of Hydrol, 143, 109-123, 1993. [10] Bogardi, I., A. Bardossy and L. Duckstein, Multicriterion network design using geostatistics, Water Resour. Res., 21(2), 199-208, 1985. [11] Bracewell, R. N., The Fourier transform and its applications, Second edition, McGraw-Hill Inc., 1986. [12] Butler, J. J. , Pumping tests in nonuniform aquifers: The radially symmetric case. J. of Hydrol., 101(1/4), 15-30, 1988. 194 Bibliography 195 [13] Butler, J. J. , A stochastic analysis of pumping tests in laterally nonuniform media. Water Resour. Res., 27(9), 2401-2414, 1991. [14] Butler, J. J. and W. Z. Liu, Pumping test in non-uniform aquifers — the linear strip case. J. of Hydrol., 128, 69-99, 1991. [15] Butler, J. J. , The role of pumping tests in site characterization: some theoretical considerations. Ground Water, 28(3), 394-402, 1990. [16] Butler, J. J. and C. D. McElwee, Variable-rate pumping tests for radially symmetric nonuniform aquifers. Water Resour. Res., 26(2), 291-306, 1990. [17] Cardwell, W. T., and R. L. Parsons, Average permeabilities of heterogeneous oil sands. Trans. AIME, v. 160, 34-42, 1945. [18] Carslaw, H. S. and J. C. Jaeger, it Conduction of heat in solids, Oxford University Press, 332-334, 1959. [19] Chirlin, G. R. and E. F. Wood, On the relationship between Kriging and state estimation. Water Resour. Res., 18(2), 432-438, 1982. [20] Cooper, H. H., J. D. Bredehoeft and S. S. Papadopulos, Response of a finite diameter well to an instantaneous charge of water, Water Resour. Res., 3(1), 263-269, 1967. [21] Cooper, H. H. and C. E. Jacob , A generalized graphical method for evaluating formation constants and summarizing well-field history. Transactions, AGU, 27(IV), August, 1946. [22] Cressie, N. A. C , Statistics of spatial data, Wiley Interscience, New York, 1991. [23] Cushman, J. H., On measurement, scale, and scaling. Water Resour. Res., 22(2), 129-134, 1986. [24] Cushman, J. H., On the unifying the concepts of scale, instrumentation, and stochas-tics in the development of multiphase transport theory. Water Resour. Res., 20(11), 1668-1676, 1984. [25] Dagan, G., Models of groundwater flow in statistically homogeneous porous forma-tions, Water Resour. Res., 15(1), 47-63, 1979. [26] Dagan, G., Stochastical modeling of groundwater flow by unconditional and condi-tional probabilities 1. Conditional simulation and the direct problem. Water Resour. Res., 18(4), 813-833, 1982. [27] Dagan, G., Analysis of flow through heterogeneous random aquifers, 2. Unsteady flow in confined formations. Water Resour. Res., 18(5), 1571-1585, 1982. Bibliography 196 [28] Dagan, G., Analysis of flow through heterogeneous random aquifers by the method of embedding matrix, 1. Steady flow. Water Resour. Res., 17(1), 107-121, 1981. 1571-1585, 1982. [29] Dagan, G., Statistical theory of groundwater flow and transport: Pore to laboratory, laboratory to formation, and formation to regional scale. Water Resour. Res. 22(9), 120S-134S, 1986. [30] Dagan, G., Flow and transport in Formations, 461pp., Springer-Verlag, 1989. [31] Delhomme, J. P., Spatial variability and uncertainty in groundwater flow parameters: A geostatistical approach. Water Resour. Res., 15(2), 269-289, 1979. [32] Demir, Z. and T. N. Narasimhan, Improved interpretation of Hvorslev tests, J. of Hydraul. Engrg., 120(4), 1994. [33] Desbarats, A. J., Spatial averaging of hydraulic conductivity under radial flow con-ditions, Math. Geol, 26(1), 1-21, 1994. [34] Desbarats, A. J. , Spatial averaging of hydraulic conductivity in three-dimensional heterogeneous porous media, Math. Geol., 24(3), 249-267, 1992a. [35] Desbarats, A. J. , Spatial averaging of transmissivity in heterogeneous fields with flow towards a well, Water Resour. Res., 28(3), 757-767, 1992b. [36] Desbarats, A. J., Numerical estimation of effective permeability in sand-shale for-mations, Water Resour. Res., 23(2), 273-268,1987. [37] Deutsch, C. V. and A. G. Journal, GSLIB, Geostatistical software library and user's guide, Oxford University Press, 1992. [38] Dykaar, B. B. and P. K. Kitanidis, Determination of the hydraulic conductivity for heterogeneous porous media using a numerical spectral approach 1. Method. Water Resour. Res., 28(4), 1155-1166, 1992. [39] Dykaar, B. B. and P. K. Kitanidis, Determination of the hydraulic conductivity for heterogeneous porous media using a numerical spectral approach 2. Results. Water Resour. Res., 28(4), 1167-1178, 1992. [40] Dong, Y.H., The groundwater quality management model of the shallow aquifer of Xinxian city, M.Sc. thesis, Dept. of Fundamental Sciences, Changchun Univ. of Earth Sci., 1994. [41] Freeze. R. A., A stochastical-conceptual analysis of one-dimensional groundwater flow in nonuniform homogeneous media. Water Resour. Res., 11(5), 725-741, 1975. Bibliography 197 [42] Freeze, R. A., J. Massmann, L. Smith, T. Sperling and B. James, Hydrogeological Decision analysis: 1. A framework. Ground Water, 28(5), 738-765, 1990. [43] Gelhar, L. W., Stochastical subsurface hydrology from theory to applications. Water Resour. Res. 22(9), 135S-145S, 1986. [44] Gelhar, L. W., Stochastical subsurface hydrology, Prentice Hall, 1994 [45] Ghilardi, P., A. K. Kai and G. Menduni, Self-similar heterogeneity in granular porous media at the representative elementary volume scale. Water Resour. Res., 29(4), 1205-1214, 1993. [46] Gomez-Hernandez J. J. and S. M. Gorelick, Effective groundwater model parame-ter values: Influence of spatial variability of hydraulic conductivity, leakance, and recharge, Water Resour. Res., 25(3), 405-419, 1989. [47] Gutjahr, A. L., Gelhar L. W., Stochastical models of spatial variability in subsurface flows: Infinite versus finite domains and stationarity. Water Resour. Res., 17(2), 337-350, 1981. [48] Gutjahr, A. L., Gelhar L. W., A. A. Bakr and J. R. MacMillan, Stochastical analysis of spatial variability in subsurface flows 2. Evaluation and application. Water Resour. Res., 14(5), 953-959, 1978. [49] Guyonnet, D, S. Mishra and J. McCord, Evaluating the volume of porous medium investigated during slug tests, Ground Water, 31(4), 1993. [50] Hassanizadeh, M., and W. G. Gray, General conservation equations for multi-phase systems, 1, Averaging procedure, Adv. Water Resour., 2, 131-144, 1979a. [51] Hassanizadeh, M., and W. G. Gray, General conservation equations for multi-phase systems, 2, Mass, momenta, energy, and entropy equations, Adv. Water Resour., 2, 191-203, 1979b. [52] Harvey, C. F., Interpreting parameter estimates obtained from slug tests in het-erogeneous aquifers, M.Sc. thesis, Applied Earth Sciences Dept., Stanford Univ., 1992. [53] Helstrom, C., Image restoration by the method of least squares, J. Opt. Soc. Am., 57, 297-303, 1967. [54] Hvorslev, M. J. , Time lag and soil permeability in ground-water observations, BULL. No.36, Waterways Exper. Sta. Corps of Engrs, U.S. Army, Vicksburg, Mississippi, 1-50. Bibliography 198 [55] Indelman, P. and G. Dagan, Upsealing of permeability of anisotropic heterogeneous formations, 1. The general framework. Water Resour. Res., 29(4), 917-923, 1993. [56] Indelman, P. and G. Dagan, Upscaling of permeability of anisotropic heterogeneous formations, 2. General structure and small perturbation analysis. Water Resour. Res., 29(4), 925-933, 1993. [57] Indelman, P, Upscaling of permeability of anisotropic heterogeneous formations, 3. Applications. Water Resour. Res., 29(4), 935-943, 1993. [58] Journel, A. G., Ch. J. Huijbregts, Mining geostatistics, Academic Press, 1978. [59] Karasaki, K., J. C. S. Long and P. A. Witherspoon, Analytical models of slug tests, Water Resour. Res., 24(1), 115-126, 1988. [60] Katz, A. J. , and A. H. Thompson, Fractal sandstones pores: Implications for con-ductivity and pore formation. Phys. Rev. Lett, 54(12), 1325-1328, 1985. [61] Keller, J. B., Darcy's law for flow in porous media and the two space method, Nonlinear partial differential equations in engineering and applied science, edited by R. L. Sternberg, A. J. Kalinowski and J. S. Papadakis, 429-443, Marcel Dekker, New York, 1980. [62] Kitanidis, P. K., Effective hydraulic conductivity for gradually varying flow. Water Resour. Res., 26(6), 1197-1208, 1990. [63] Marie, C.-M., Ecoulements monophasiques en milieu poreux, Rev. Inst. Fr. Petrol., 22, 1471-1509, 1967. [64] Marsily, G. de, Quantitative hydrogeology: Groundwater hydrology for engineers, Academic Press, Orlando, FL, 1986. [65] Matheron, G., Les Variables Regionalisees et Leur Estimation, Masson, Paris, 1965. [66] Matheron, G., Elements pour une Therie des Milieux Poreux, Masson et Cie, Paris, 1967. [67] McElwee, C. D., Sensitivity analysis of groundwater models, Advances of Transport Phenomena in Porous Media, J. Bear et al. (eds.), NATO adv. Study Inst. Ser., Ser.E, 128, 751-817, 1987. [68] McElwee, C. D., M. A. Yukler, Sensitivity of groundwater models with respect to variations in transmissivity and storage, Water Resour. Res., 14(3), 451-459, 1978. Bibliography 199 [69] Moench, A. and A. Ogata, Analysis of constant discharge wells by numerical inver-sion of Laplace transform solutions. Groundwater Hydraulics. Am. Geophys. Union Water Resour. Monogr.9 , J. Rosenshein et al. (eds.), American Geophysical Union, Washington. DC., 146-170, 1984. [70] Mohanty, N., Random signal estimation and identification, Van Nosttrand Reinhold Company Inc., 1986. [71] Murtagh, B. A. and M. A. Saunders, MINOS, A modular in-core nonlinear opti-mization system, Stanford University, 1983. [72] Naff. R. L., Radial flow in heterogeneous porous media: An analysis of specific discharge. Water Resour. Res., 27(3), 307-316, 1991. [73] Newland, D. E. , An introduction to random vibrations, spectral & wavelet analy-sis(Third edition), Longman Scientific Technical, 1993. [74] Oliver, D. S., The influence of nonuniform transmissivity and storativity on draw-down, Water Resour. Res., 29(1), 169-178, 1993. [75] Oliver, D. S., Estimation of radial permeability distribution from well-test data, SPEFE, 290-296, December, 1992. [76] Oliver, D. S., The averaging process in permeability estimation from well-test data, SPEFE, 319-324, September, 1990), [77] Oppenheim, A. V-, A. S. Willsky, and I. T. Young, Signals and systems, Prentice-Hall, Inc., 1983. [78] Papadopulos, S. S., On the analysis of 'slug test' data, Water Resour. Res., 9(4), 1087-1089, 1973. [79] Park, J. , C. R. Lindberg and F.L.Vernon, Multitaper spectral analysis of high fre-quency seismograms, J. geoph. Res., 92, 12675-12684, 1987. [80] Peaceman, D. W., interpretation of well-block pressures in numerical reservoir sim-ulation, Society of Petroleum Engineers of AIME, 6, 183-194, 1978. [81] Priestley, M. B., Spectral analysis and time series, Vols. I h II, Academic Press, N.Y., 1981. [82] Rubin, Y. and J. J. Gomez-Hernandez, A stochastical approach to the problem of upscaling of conductivity in disordered media: Theory and unconditional numerical simulations. Water Resour. Res., 26(4), 691-701, 1990. Bibliography 200 [83] Sageev, A., Slug test analysis, Water Resour. Res., 22(8), 1323-1333, 1986. [84] Sanchez-Vila, X, J. P. Girardi and J. Carrera, A synthesis of approaches to upscaling of hydraulic conductivities, Water Resour. Res., 31(4), 867-862, 1995. [85] Schwydler, M. I., Flow in heterogeneous media(in Russian), 185pp., Gaza, 1962. [86] Schwydler, M. I., Sur la precision de la precision du debit des piuts, ize, Akad. Nauk SSSR Mekh. Mashinostr., no. 5, 148-150,1963. [87] Smith, L. and R. A. Freeze, Stochastical analysis of steady state groundwater flow in a bounded domain, 2. Two-dimensional simulations. Water Resour. Res., 15(6), 1543-1559, 1979. [88] Sondhi, M. M., Image restoration: the removal of spatially invariant degradations, Proc. IEEE, 60, 842-853, 1972. [89] Stehfest, H., Numerical inversion of Laplace transforms, Commun. ACM, 13(1), 47-49, 1970. [90] Sun, N. Z. and W. W-G. Yeh, Identification of parameter structure in groundwater inverse problem, Water Resour. Res., 21(6), 869-883, 1985. [91] Theis, C. V., The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground-water storage. Am. Geophys. Union Trans., 16:519-524, 1935. [92] Thomson, D. J., Spectrum estimation and harmonic analysis, Proc. IEEE, 70, 1055-1096, 1982. [93] Vanmarcke, E. , Random Fields; Analysis and Synthesis, MIT Press, Cambridge, Mass., 1983. [94] Warren, J. E. , F. F. Skiba and H. S. Price, An evaluation of the significance of permeability measurements. SPEJ, August, 1961. [95] Wang, B., Groundwater pollution and mathematical modeling in the shallow aquifer of the chemical-industry district of Taiyuan city, PRC. M.Sc. thesis, Chinese Acad, of Geol. Sci., 1992. [96] Yeh, W. W-G., Review of parameter identification procedures in groundwater hy-drology: The inverse problem, Water Resour. Res., 22(2), 95-108, 1986. Appendix A E F F E C T I V E A N D E Q U I V A L E N T P A R A M E T E R S We have illustrated in eq.(1.5) that the real world has already been mapped to an en-gineering framework after the measurement filtering HP and estimation filtering HS. However, HS filtering is a probabilistic or stochastic process, and it generates infinite number of realizations in which each realization has equal significance. These realiza-tions may require exorbitant computation load when they are utilized in the simulation of engineering problems(e.g., Monto Carlo simulation), even though they may well suit the scale-requirement of the decision variable. The consequent compromise is the use of effective parameters, which are defined under ergodicity as the representative values of the mean behavior through an ensemble of realizations [Sanchez- Vila et al., 1995]. Equivalent parameters are defined as spatial averages obtained on a single realiza-tion [Sanchez- Vila et ai, 1995], not necessarily requiring ergodicity or stationarity. The equivalent parameters are associated with certain initial and boundary states. Conse-quently, equivalent parameter is deterministic. If the equivalent parameter scale equals to the size of a grid block, we call it a block equivalent parameter, or shortly, a block parameter. Apparently, block parameters are dependent upon the block geometry and block boundary states, but this dependence can be lifted under certain conditions, e.g., increasing the block size(reducing the boundary effects) and imposing stationarity on the upper-block scale. Indeed, under stationarity and ergodicity, when the block size is enlarged to infinity, or more strictly, the block scale is large enough relative to global het-erogeneity scale, the block parameter approaches the effective parameter. Consequently, 201 Appendix A. EFFECTIVE AND EQUIVALENT PARAMETERS 202 block parameter a generalization of the effective parameter. A.I Effective parameter: a limiting case of upscaling The engineering importance of the real world heterogeneity has motivated the study effective parameters, of which conductivity receives the most concern. Among the pi-oneering work on several disciplines, Warren and Price [1961], using a Monte Carlo method, studied the flow in heterogeneous media in petroleum engineering; Matheron [1967] relates the average permeability to the harmonic, geometric and arithmetic mean of the local permeabilities, in the geostatistics context. The Monto Carlo simulation of Freeze [1975], followed by Smith and Freeze [1979] led to the explosion of heterogeneity study from various aspects in water resources research. Following the pioneering works, enormous efforts have been devoted to the study of effective conductivity in heterogeneous media(hereafter Kej denotes effective conductiv-ity). With both numerical and analytical vehicles, one is able to relate the statistics of Kef to those of small scale conductivities assumed to be known. Analytical tools include the versatile first-order perturbation method [Schwydler, 1962; Matheron, 1967; Bakr et al, 1978; Gutjahr et al., 1978; Gutjahr and Gelhar, 1981], self-consistent method [Dagan, 1979; 1989], and the method of volume-averaging and spatial moments [Kitanidis, 1990]. In order to account for more complicated flow and porous media and initial and boundary effects, another group of authors used semi-analytical method [Desbarats, 1992a; 1992b; 1994], and Monto Carlo method [Freeze, 1975; Warren and Price, 1961; Desbarats, 1987; Gomez-Hernandez and Gorelick, 1989]. The review shall not be exhaustive but comprise the major contributions. Other more detailed summary can be found in Gelhar [1994], Gomez-Hernandez and Gorelick [1989], Sanchez-Vila et al. [1995], and Desbarats [1987, 1992a]. Appendix A. EFFECTIVE AND EQUIVALENT PARAMETERS 203 Under the condition of uniform steady mean flow in infinite isotropic media of small variance, the effective conductivity Kej equals to K E F = KGex7[<7(\-l)l~ ( A - L 1 3 5 ) where KG is the geometric mean and n = 1, 2, 3 are flow dimensions [Gutjahr et ai, 1978; Bakr et al, 1978; Gutjahr and Gelhar, 1981; Dagan, 1979; 1989]. The 2D solution was extended to any CTJ^ by Gelhar [1986]. The 3D solution is the first-order approximation of the Kej first postulated by Matheron [1967]. ~ In transient flow Freeze [1975] states Kef is time-dependent. In the slowly varying flow excited by a instantaneous point 'slug injection' in conductivity field of small C T k ( F = In if), Dagan [1982] found that Kej varies from arithmetic mean at initial time to steady state value during relaxation time(which must elapse after some change in the system so that Kef can be defined). The flow condition is similar to that of a slug test of Cooper et al. [1967], except that K and Ss at the injection location are of that of aquifer as usual, whereas at the wellbore of a slug test K = oo, Ss = l . 1 Kitanidis [1990] studied the effective conductivity in the flow condition similar to Dagan [1982] but in periodic conductivity fields. Kitanidis [1990] derived Kef by the generalization of block conductivity Key? He represented Kev in a matrix form(eq.56), ( A i ) e V = - 2 ^ h h j [di • A9j + dj • A9i] dx+ D~j, (A.I.136) where D is diffusivity(K/5), overbar denotes volumetric average, / is the side lengths of rectangular parallelepipeds, and g is a function of local spatial coordinates. (Kij)eV can be obtained by solving the matrix if Sa is constant. In a isotropic media, above equation 1 Dimensions of K and S are [LT - 1] and [L - 1] respectively, they will be omitted in the following discussion. 2Kev will be discussed in next section. Appendix A. EFFECTIVE AND EQUIVALENT PARAMETERS 204 reduces to where, K = ^ K(x)dx (A.I.138) f>ij is the Kronecker delta, <5,j = 1 if i = j and Sij = 0 otherwise. For perfectly stratified media, above further entails (tfn)ev = (K22)eV = K = KA, (A.I.139) (•^33)ev = 1 / du k J K(u) (A.I.140) K(u) where KA denotes arithmetic average. For the media with small variations, the volume-averaging effective conductivity can be expressed in a general equation. {K^ = -exp ( F ) £ ^ |*(fc)|2 + exp ( F ) [l + 4/2] 6i};(A.I.141) where ^(A:)|2 is the power spectrum of the periodic function, and k2 = k\ + k\ + k\. If V —• co, Kev —• Kej, by replacing the summation of a integral, above equation reduces to • (Kij)eJ = -exp ( F ) J ^S(k)dk + exp ( F ) [l 4- a2Y/2] ^ (A.I.142) where S(k) is the PSD of Y. This is identical to the result of Gutjahr et al. [1978] which followed a ensemble — averaging approach, whereas here follows a purely deterministic volume—averaging approach. This essentially means the result for a periodic medium by taking a limit(the period tending to infinity) is applicable to a random medium. (K(j)ej can be reduced to those of Gutjahr et al. [1978], Bakr et al. [1978], and Gutjahr and Gelhar [1981], also confirms the 3D effective conductivity hypothesis of Matheron [1967], but involving no assumption other than isotropic conductivity variability. Appendix A. EFFECTIVE AND EQUIVALENT PARAMETERS 205 Kitanidis [1990] showed the assumptions of steady flow [Gutjahr et al., 1978; Bakr et al., [1978]; Gutjahr and Gelhar, 1981] and unsteady flow [Dagan, 1982]) are not essential in the definition or determination of effective parameters. Instead, the essential assump-tion is that of 'gradually varying' flow. That is the scale of fluctuations of the head must be large in comparison to the typical length scale of fluctuations of local conductivity so that all values of conductivity are sampled(eq.47). Kitanidis [1990] also found that shortly after the slug injection the net spreading rate of slug mound depends on the value of conductivity in the neighborhood of the injection and varies with time, and this is when the 'gradually varying' condition is not satisfied, hence, this is when effective conductivity is not meaningful. This spreading continues up to the relaxation time and converges to a constant, and this time can be approximated by P SalK, where 7 is the scale of conductivity fluctuations, K is a typical value of conductivity. These observations are in general agreement of Dagan [1982]. In his study of 2D radial flow in an annular region drained by a pumping well with fixed head drawdown at (r = rw) and zero drawdown at some 'large distance' r = r m , Matheron [1967] found that in the unconditional case, Kef —• KH for rm/\ —• oo; and Kef —> KA for rw/X —> 0; while in the conditional case — with K = Kw at the well, Kej —* KH for r m / A —• oo; and Kef —* Kw for rw/\ —* 0, where A is conductivity correlation length, KA and KH are the arithmetic and homonic mean, respectively. Gomez-Hernandez and Gorelick [1989] performed Monto Carlo simulations under a quite complicated conditions: unconfined aquifer in bounded domain with a variety of boundary conditions, and distributed recharges and recharge from river, and having pumping wells, they found KH < Kef < K Q . They used the Kej defined in p norm to fit the Kef calculated from MC simulation, Kp = (Kp)1/p (A.I.143) Appendix A. EFFECTIVE AND EQUIVALENT PARAMETERS 206 The numerically determined p values are confirmed in the analytical derivation of Ababou and Wood [1990], directly relating to the statics of Y(= In K), KG exp (p 0Y/2) ~P'= KG (KG/K„Y, p = 0, (A.I.144) [KG(KA/KGy, p = - l Ababou and Wood [1990] also expanded the 2D Kef results of Matheron [1967] by com-bining the perturbation solutions from Matheron [1967, section 25] and Schwydler [1963]. They found the effective conductivity for a Gaussian-shaped correlation function as K ef KGexp 1-2- _(ln(rm/rw)Y = r^exp |y^2 [l + ( ln^ /A) ) 2 ] ! With the assumption rw «C A <C rm, it yields -11 (A.I.145) (A.I.146) Kef < KG, rm > rG Kef > KG, rm<rG (A.I.147) [ Kef = KG, rm = rG where rc is a critical depression cone radius called "geometric'1 radius. Besides above findings, they stated: (i) For fixed A, Kef —• Kn, for r m —• oo; Kef decreases and becomes closer to KH as rw increases but remains smaller than A; (ii) Kef —* K A , for rw —¥ 0; Kef increases as rm decreases but remains larger than A; and (iii) For fixed rw and rm, the Kef increases as A increases within [»*„,, rm]. This somehow explains the Kfj < Kef < KG of Gomez-Hernandez and Gorelick [1989]. Note that Dagan [1979] also found this Kn < Kef < KG relationship in unbounded isotropic fields of small variance, but Kef will exceed KG for larger variances. Using first order perturbation analysis, Naff [1991] studied the 3D steady mean-radial-flow to a pumping well sited in infinite medium. Naff [1991] found: (i) At distance Appendix A. EFFECTIVE AND EQUTVALENT PARAMETERS 207 p = r/A > 2 or 3, if e/ remains constant, depending on the anisotropy p = A/A3; (ii) At the wellbore, Kef = K A , regardless of p; (iii) At large distance, Kej varies from A#(2D flow) to KA(SD flow in a well-stratified medium). In general the results agree closely with the asymptotes of Matheron [1967]. Interestingly, Naff [1991] stated the variance of specific discharge^/) is approximately proportional to r - 2 . 3 A . I I Equivalent parameter: a general case of upscaling As characterized by Sanchez-Vila et. al [1995], there are four major types of block up-scaling methods reported in literature: (i) Darcian definition; (ii) Method of Rubin and Gomez-Hernandez [1990]; (iii) Power averaging method of Desbarats [1992a]; and (iv) Method of Indelman and Dagan [1993a, b] and Indelman [1993]. Sanchez-Vila et. al [1995] state that all the methods seem to provide striking agreements in statistic terms. Here we only discuss method (i) and (iii), and one more method - weighted averaging method [Dykaar and Kitanidis, 1992a, b] The upscaling by the Darcian approach directly applies to the equation of fluid flow in porous media at small scale V. [ * ( * )V* ( * ,0 ] = S ^ J T - , (A.II.148) upscaling to a large-scale equation, where the small-scale effects are lumped into the effective parameter Kev- The proportionality of a upscaled Darcy's law equals **• - W (A-IU49) and the upscaled flow equation is V-[KeV Vt f * , * ) ] = S ^ ^ . (A.II.150) 3 T h i s may be somehow explained by the r " 2 in Kfv = & fv w = Sv r T ^ F ° f D e s b a r a t s [1994], where V is the drainage volume of pumping test and *" is an empirical exponent. By applying Darcy's law, qef(r) = Kej(r) Jej(r) where J is gradient, the r ~ 2 reaches qef. Appendix A. EFFECTIVE AND EQUIVALENT PARAMETERS 208 The upscaling should meet the scale requirement [Dagan,19S6; Gelhar, 1986], / < L < D (A.II.151) where / is the small-scale, L is the scale of the averaging volume, and D is the global heterogeneity scale. Desbarats [1992a] empirically defines the effective conductivity as the "power average" of random function K\x) over a volume V KeV = (1 JvK(x)wdVy (A.II.152) and analyzes the exponent w as that best fits numeric simulations. Dykaar and Kitanidis [1992a, b] defines Kev in as a weighted spatial average in a volume V, KeV(xc) = + (A.II.153) v where xc is the centroid of V and W(£) a weight function matrix obtained by solving flow equations. This is similar to the definition of the measurement filter of eq.(1.2) where the filter is equivalent to the weight function. The idea of weighted average was described by a number of authors [Matheron, 1965; Marie, 1967; Anderson and Jackson, 1967] (cf. the discussion in Baveye and Sposito [1984]). Dykaar and Kitanidis [1992a, b] found that in a 2D flow in conductivity fields of a Gaussian model with aY — 2, 4 and 5 < L/l < 80, a V of 80A(integral scale) is needed for Kev to reach a relative constant value, while in 3D a V of 30A is required for Kev to reach its asymptotic value. The heterogeneities of scale smaller than 2A(for 2D) and 1.3A(3D) do not significantly contribute to Key-Appendix B F I N I T E D I F F E R E N C E E Q U A T I O N S O F S L U G T E S T B . I 2D F D M in rectangular coordinates Forward finite difference approximation of equation (2.16), Figure B.I.80, dxx • i,j+l W.J \ -8-Figure B.I.80: Block centered finite difference in Cartesian coordinates. m 1 hm — hm. hm — h j- [ K i j H ( ^ M ± I _ A i ) _ K ( K ± ^ M ) ] uxx Uxp d i m 1 hm • — hm hm — hm + d 1 *+*•''1 5 — } ~ K i - y ( — 1 — ' ] "yy u y p "ym 209 Appendix B. FINITE DIFFERENCE EQUATIONS OF SLUG TEST 210 AT™. _ A ™ . - 1 = S. *'J , , J (B.I.154) Multiplying by dxxdyy, to get the mass balance equation at block (i, j), A™. _ A m . A"» _ l m 3 ) - Ki,i-\\ 3 )J A m . — A m . A"* _ A m uyp " aym fl™. — / j ? 1 . - 1 = dxxdyySa-^-^— (B.I.155) Rearrange the above, to get, where, K { i . l K;-L K;,l- K; 1 • A A C Lv J + ~3 + \—j r —; ja** + —Jft.j l*xp 8 i m "yp "ym u xp ttim = dxxdyySs, m _ i A* 6?xp — rfx(j) + dx(j - 1) 2 dx(j) + + 1) 2 dXm ~T" ^xp o dym = 2 <fy(i) + - 1) 2 dyp = <*y(«) + <*y(« +1) 2 (B.I.156) Appendix B. FINITE DIFFERENCE EQUATIONS OF SLUG TEST 211 K. (dyy + dyp)Ki+ijKij dypKij + dyyKi+lj I i = ' dymKij + dyyKi-lj K, *'3+2 dxpK~itj + dxxK~itj+i _ (dxx + dxm)Kij-\Kjj 2 dxmK{tj -f- dxxKij—\ Kij-h ~ J is i J is — (B.I.157) Finite differencing wellbore equation (2.17), yields the discretized mass balance equation at Wellbore, (ij) = (iwellborejwellbore), hm. — hm. hm — hm Syl A.J+|( ) - K i J - \ \ J £ )J xp xm hm . — hm. hm — hm + dr \K. i ( , + 1 , J - K i ( , J l-1'J)] At = "J A / " (B.I.158) Equation (B.I.158) is the same as equation (B.1.155) except Ss = 1 in the right hand side term. This physically means that wellbore is filled with water. Equation (B.I.158) is the same as equation (B.I. 155) except Ss = 1 in the right hand side term. This physically means that wellbore is filled with water. Rearrange the above, get, Kw.^_ Kw. , Kw. . Kw1 . dw dw [ ( - ^ i + _ ^ = I ) d £ , + ( ^ ± M + + -^*L]kT. ux p "xm Uyp "ym L X L dw dw ( yy Tsw u m ( yy vw \Lm xp Appendix B. FINITE DIFFERENCE EQUATIONS OF SLUG TEST 212 A V = Ax Ay AH Ss = Ax Ay AH Ss = 1 Wellbore specific stoarge equals to 1. dw dw ("•xx jy-w u m ( "xx ryw u r n VP ym dw dw _ uxxayy • m-1 " At ^ Equation (B.I .159) is the same as equation (B.I.156) except Sa = 1 and (B.I.159) <c = dx(j - 1) 0 d™ = xp dx(j + 1) 2 dwTT = XX rfx(j) < m = dy(i - 1) 2 rfy(i + 1) 2 <C = dy{i) (B.I.160) Appendix B. FINITE DIFFERENCE EQUATIONS OF SLUG TEST 213 (B.I.161) where, w indexes coefficients associated with wellbore (i,j). To keep symmetry, need modify equations at nodes to the east, south, west and north of wellbore East node (ie, je) = (i,j + 1) d. dx(je) xm " I I — "xm i 2 dfji = K(ie,je) South node (ia,ja) = (i — l,j) dyp dy(ia) 2 1 2 d i a + i j t = K(ia,ja) j I ( "yp West node (iw, jw) = (i,j - 1) dXp dx(jw) 2 « x x — «xp T 2 d f „ j . + i = K(iw,jw) (B.I.162) (B.I.163) (B.I.164) Appendix B. FINITE DIFFERENCE EQUATIONS OF SLUG TEST 214 North node (i n,Ji n ) = (« + W ) <fy(en) 2 (B.I.165) Coefficient definitions from (B.1.160) to (B.I.165) reflect attempts to approximate the inner boundary which encloses a non-aquifer area filled with water. Comparing the corresponding coefficient terms in (B.I.156) and (B.I.159), it immediately shows that if setting Ss in (B.I.156) equal to 1 and conductivity at block several magnitude greater than that of surrounding blocks, (B.I.156) collapses to (B.1.159). This is a much easier way to implement inner boundary. Another big advantage of this technique is that it makes it feasible to assign any number of blocks as wellbore so as to get a finer wellbore discretization, thus reducing numerical error, see Section 2.2.1. B.II 2D F D M in radial coordinates Mesh-centered finite difference equation, Figure B.11.82, rl} A0 ) ) = S, m-1 (B.II.166) s Multiplying by r l . \ { . to obtain symmetry, Appendix B. FINITE DIFFERENCE EQUATIONS OF SLUG TEST 215 Figure B.11.82: Mesh centered finite difference in radial coordinates. Ki+i ,• r . - . i K{_i ,• r . _ i A / S rl Al-. r . _ i_ i 7 ' r/, A 0 2 K _ i , r{_i 2" ' 2 Um _ ' 2 ' J ' 2 i m ' » 2 LT71 A " - " ' r / , A 02 -I 5 S W, A/,- j A* 1 7 m— (B.II.167) where, 2 Kjj+iKjj 2 Kjj-iKjj K: (A/ , + 1 + Ali)Ki+1jKij Ali+1Kitj + AUKi+ij » 2 ' J (A/,- + Alj-^KjjKj-ij AliKi-hj + Ali-iKij Appendix B. FINITE DIFFERENCE EQUATIONS OF SLUG TEST 216 r(i) + r(i — 1) ' 2 2 A/,-A / , + 1 r(Q + r(» + 1) 2 r( i + 1) — r(i) = r ( i ) - r ( » - l ) r ( i + 1) - r(* - 1) r/,- = r , _ i + — 2 2 2 r(t + 2) - r ( i ) 2 M i = r(«) - r ( i - 2) ( B Figure B.11.83: Finite difference of the inner boundary. Appendix B. FINITE DIFFERENCE EQUATIONS OF SLUG TEST 217 r c A0 Ki hm • hm- - hm-CA0 — rB rcA0 Multiplying by r;,, rearrange, coupling to domain flow equation, n K*h<' hm _ r ° T b KiJ+h hm _ r ° V b K i d - \ h m / . + 4 r 2 A 0 2 r 2 A 0 2 feT1,"1 (B.II.170) r b r c tra - 1 2A* i J where, ra = n = r(2) + r(l) 2 r(2) - r(l) Kij+i,Kij-i,Ki+ij have same definition as (B.II) except A/,- = which means that wellbore is a non-aquifer area. B.III ID F D M in radial coordinates The finite difference equation, Appendix B. FINITE DIFFERENCE EQUATIONS OF SLUG TEST 218 Rearrange the above, - - - (B.III.172) K At where, = r ( t + 1) - r ( » ) = r ( t ) - r ( t - l ) r( i — 1) + r(i) r,- i = — - — i - 5 2 = r ( Q + r ( » + 1) «+| 2 A/,- r ( t + 1) - r(i - 1) 2 r/ ,- = r . - . i + ^ i (B.III.173) Finite differencing wellbore equation, hm — hm hm — hm~l 27T rc K h>+] = x r ? ^ - ^ -c Ai Rearrange, get where has same definition as (B.III.173). (B.III.174) ( 2 K + K - 2 K h?+1 = h?-1 (B.III.175)
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The spatial filtering from core-scale conductivity...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The spatial filtering from core-scale conductivity to the conductivity measured by a slug test Wang, Bin 1995
pdf
Page Metadata
Item Metadata
Title | The spatial filtering from core-scale conductivity to the conductivity measured by a slug test |
Creator |
Wang, Bin |
Date Issued | 1995 |
Description | In the thesis, I investigate the measurement process of the slug test, an important and efficient engineering technique employed to probe aquifer properties. The objective is to determine the scale at which a slug test 'measures' hydraulic conductivity. To determine the scale, I conceptualize slug test measurement process in a system and filtering framework. I quantify, through numeric experiments, the spatial filtering effect of the slug test upon the core-scale hydraulic conductivity field. I examine the role of specific storage on this filtering effect. I develop two approaches to identify the system filter that epitomizes the filtering system — measurement process: spectral analysis with Wiener filtering and numerical perturbation. In heterogeneous media, system linearity is analyzed and the inverse estimation of the system output is evaluated. With the spectral approach, I optimally estimate the system filter in the presence of measurement noise. In homogeneous media where the spectral approach cannot be used, I apply the numerical perturbation approach. The spectral analysis and numerical perturbation are two conceptually different but complementary approaches. I analyze the nonparametric form of the system filter in heterogeneous media of various spatial variability and specific storage values. I determine the empirical parametric expression of equivalent filter width versus specific storage in homogeneous media. It is found that: 1. In both heterogeneous and homogeneous media, the equivalent filter width increases linearly with respect to log [mathematical formula which cannot be reproduced here; see pdf] 2. In homogeneous media, the filter amplitude decreases away from the wellbore by an approximate [mathematical formula which cannot be reproduced here; see pdf]. 3. A preliminary analysis of the role of heterogeneity characteristics upon the slugtest filtering was undertaken. The slug-test filtering was not strongly influenced by heterogeneity characteristics for log-conductivity variances less than 0.7. At variances less than 0.3, anisotropy in the conductivity field caused anisotropy in the slug-test filtering. These effects were not quantified in this thesis. In the end, I discuss the possible engineering and theoretical applications of the spatial filtering approach and make recommendations for future work. |
Extent | 14237754 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-02-09 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052911 |
URI | http://hdl.handle.net/2429/4352 |
Degree |
Master of Applied Science - MASc |
Program |
Geological Sciences |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1995-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1995-0528.pdf [ 13.58MB ]
- Metadata
- JSON: 831-1.0052911.json
- JSON-LD: 831-1.0052911-ld.json
- RDF/XML (Pretty): 831-1.0052911-rdf.xml
- RDF/JSON: 831-1.0052911-rdf.json
- Turtle: 831-1.0052911-turtle.txt
- N-Triples: 831-1.0052911-rdf-ntriples.txt
- Original Record: 831-1.0052911-source.json
- Full Text
- 831-1.0052911-fulltext.txt
- Citation
- 831-1.0052911.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052911/manifest