- © 2003 Society of Exploration Geophysicists. All rights reserved.
We describe the implementation of a versatile method for interpreting gravity and magnetic data in terms of 3D structures. The algorithm combines a number of features that have proven useful in other algorithms. To accommodate structures of arbitrary geometry, we define the subsurface using a large number of prisms, with the depths to the tops and bottoms as unknowns to be determined by optimization. Included in the optimization process are the three components of the magnetization vector and the density contrast, which is assumed to be a continuous function with depth. We use polynomial variations of the density contrast to simulate the natural increase of rock density with depth in deep sedimentary basins. The algorithm minimizes the quadratic norm of residuals combined with a regularization term. This term controls the roughness of the upper and lower topographies defined by the prisms. This results in simple shapes by penalizing the norms of the first and second horizontal derivatives of the prism depths and bottoms. Finally, with the use of quadratic programming, it is a simple matter to include a priori information about the model in the form of equality or inequality constraints. The method is first tested using a hypothetical model, and then it is used to estimate the geometry of the Ensenada basin by means of joint inversion of land and offshore gravity and land, offshore, and airborne magnetic data. The inversion helps constrain the structure of the basin and helps extend the interpretation of known surface faults to the offshore.
In geophysics, the gravity method seems to have the most variety of techniques devised for data interpretation, in particular what we now call inverse methods. The techniques range from the old but clever half-width rules for depth determinations of spheres and cylinders (e.g., Nettleton, 1940) to rigorous theories for establishing true bounds on density variations (e.g., Parker, 1974). Somewhere in between are the popular data fitting algorithms that automatically provide models whose responses reproduce the data to reasonable levels. These fitting algorithms differ from each other in several ways. One of them is the particular way that density variations are parameterized. In this respect, we can distinguish two groups: in one, the inverse problem is linear; in the other, it is nonlinear. In the first group, the recovery of the density distribution is linear because this is how it appears in Newton's gravitational law, where density is assumed as an arbitrary function of space (e.g., Green, 1975). In the second group, the model consists of one or more bodies of uniform densities and unknown geometry, in which case the combined recovery of densities and shapes leads to a nonlinear problem (e.g., Cordell and Henderson, 1968). The payoff for opting for the nonlinear alternative is that the nonuniqueness of the inverse problem is less severe than in the linear option. Silva et al. (2002) and Barbosa et al. (2002) recently reviewed this and other aspects of uniqueness in relation to gravity inversion. The reader is referred to these papers for an extensive bibliography on the subject.
Given the great variety of available approaches for inverting potential field data, it would seem that when faced with the problem of interpreting a given set of gravity and magnetic data, all we need to do is simply select the appropriate approach from the pool of existing methods. Unfortunately we are not there yet, at least not for 3D situations. This explains the large number of papers on the subject and the fact that it is still an active area of research. Consider, for instance, the problem of jointly interpreting a set of data that includes land, airborne, and marine gravity and magnetic measurements taken at different altitudes, with a 3D model of a partly offshore, partly inland sedimentary basin, whose upper and lower boundaries are arbitrary functions of position, and whose mass density contrasts vary with depth to accommodate possible density variations due to natural loading. We could not find in the literature a method to handle this situation. Some methods can handle both gravity and magnetic data but do not handle variable density (e.g., Zeyen and Pous, 1993). Some others deal with variable density but cannot handle data for different altitudes (e.g., García-Abdeslem, 1992), or they limit the model to the lower boundary only (e.g., Chai and Hinze, 1988).
The problem calls for a type of algorithm that combines features of both linear and nonlinear inversion. On one side, we need variable density limited to the vertical dimension; on the other, we need a structure with very definite boundaries. Both features can be accommodated by assembling a set of vertical prisms for defining the upper and lower boundaries of structures and by assuming a simple quadratic law for the vertical variation of density. The combination of these two features has been exploited by some of the authors cited above, but in a simplified form that prevents the inclusion of data taken at different levels. The generalization of existing formulas to handle multilevel data remained to be done and applied to inversion, if for no other reason than to avoid tedious data merging and re-leveling.
Another weakness of existing algorithms concerns the proper stabilization of the geometry of the model, which worsens when going from one to two bounding surfaces. Applying regularization through the second spatial derivative of the boundaries solves this problem effectively. This technique, popular in electromagnetic inversion, is adapted here from derivatives of physical properties (e.g., DeGroot-Hedlin and Constable, 1990) to derivatives of geometry. Still another aspect that existing methods have difficulty to deal with is the application of equality and inequality constraints, a tool extremely useful for the inclusion of a priori information. This we handle here with the use of quadratic programming.
First, we present an analytical expression for the gravity effect at any altitude from a prism whose density varies with depth in a quadratic fashion. This expression is a simple but nontrivial extension of existing formulas and one of special utility in a wide variety of applications. Without it, it is not possible to invert data taken at different altitudes, as the case may be of land and aerial data. Neither can the lower topography of a basin be determined when the data is taken over irregular upper topographies. Even with data taken at the same level, as the case may be with marine gravity, we require the general expression for modeling offshore basins with irregular sea floors. The magnetic effect of a prism is also presented, followed by the development of analytical partial derivatives for both gravity and magnetic fields with respect to all the relevant parameters. Next, we provide a brief description of the optimization technique, and an application of the inverse algorithm to the recovery of a relatively simple model. Emphasis is placed on the capabilities of quadratic programming for placing inequality constraints, a most important feature for managing inherent ambiguities. Finally, we present the inversion of the combined set of data that we could not interpret with existing algorithms. The application determined the geometry of the offshore Ensenada Basin and helped to identify previously unknown faults.
GRAVITY AND MAGNETIC RESPONSE OF A 3D PRISM
Newton's law of gravitation for the vertical component of gravity gz at a point (x0, y0, z0) due to a rectangular prism can be written as (1) where ρ is the gravitational constant, and r is the distance between (x0, y0, z0) and (x, y, z) as represented in Figure 1. The limits in the integral signs define the prism dimensions; ρ represents the mass density function, which in the present case we assume as (2) where a, b, and c represent constants. We assume that the density does not change laterally. Substituting this density law into equation (1) and performing the corresponding integrations, we arrive at (3) The limits Δ xi = xi ∼ x0, Δ yi = yi ∼ y0, and Δ zi = zi ∼ z0 are to be evaluated at i = 1, 2 as defined in Figure 1. The depth variable z is evaluated at the top, z1 = ht, and bottom z2 = hb of the prism.
In the derivation of equation (3), we followed the steps of Bhaskara-Rao et al. (1990), who obtained the corresponding expression for z0 = 0. Notice that for z0 = 0, the polynomial factors a + bz0 + cz02 and b + 2z0 c, which stand in front of the first two bracketed terms in equation (3), cancel any information about the altitude of the measurement. The resulting expression is the one derived by Bhaskara-Rao et al. (1990), and corresponds to the special case when the measurements are taken exactly over the prism. It is not possible to simply change variables in their result to obtain equation (3); the integration needs to be redone considering expression (2), which defines the density variations within the prism from arbitrary altitudes. The required integrals are very similar to the ones for the case of z0 = 0, and we will not present the full development here.
In the case of the magnetic response of a 3D prism, we use the total component of the magnetic field, which can be expressed as (4) where J = (Jx, Jy, Jz) represents the magnetization vector, which up to this point may vary with position; k = 1 in emu units. The rest of the variables are as in equation (1).
equations (3) and (5) are the basic expressions for computing the effect of a prism at any altitude: the first for the vertical component of gravity and the second for the total component of the magnetic field. The response of complex bodies made up of elementary prisms as shown in Figure 2 are computed simply by adding the individual responses. equations (3) and (5) are also the starting points for developing the partial derivatives with respect to the unknown quantities. The resulting expressions are developed in the Appendix.
THE FITTING ALGORITHM
We need to find a model whose gravity and magnetic responses reproduce to a reasonable level a given set of observed and magnetic data. The model is characterized by the three coefficients a, b, and c of the density law, the three components of the magnetization vector, and the set of top (ht) and bottom (hb) depths of the array of prisms depicted in Figure 2. It is clear from equations (3) and (5) that the three coefficients and the components of the magnetization vector are linearly related to their corresponding responses. It is also clear that the relationship of the responses with the top and bottom depths is of a nonlinear type. It then follows that any fitting algorithm for the simultaneous estimation of all parameters must necessarily be nonlinear. Still, as is usually done in these cases, it is convenient to pose the problems as a linear one. To do this, we begin by considering a standard Taylor expansion of gz from equation (3) around a reference model. Neglecting second-order terms, we have that gzj, the jth gravity measurement, can be written as (6)
The vector of unknown variables X = (X1,…,X2n + 3) includes the three coefficients a,b,c of the density function, and the n top (hit,i = 1,n) and the n bottom (hib,i = 1,n) depths of the n prisms. The superscript “0” in the vector X0 = (X10,…,X2n + 30) indicates that the variables are to be evaluated at the reference model. The term gzj(X0) is the jth gravity response of the reference model at the same location and altitude of the corresponding gravity measurement. The derivatives of gzj(X0) with respect to the relevant parameters are derived analytically from equation (3) in the Appendix.
For a given set of measurements, equation (6) represents a system of linear equations whose solution is (X ∼ X0). Even though X, the vector of unknowns, can be computed from this difference, it is more convenient to solve directly for X, because it is easier to impose constraints on the unknown parameters in this way than through their difference. With this in mind, equation (6) can be rearranged as (7)
For a number mg of gravity measurements, we have a system of mg equations with 2n + 3 unknowns. This system can be written as (8)
In the case of magnetic data, we can arrive to an expression similar to equation (7), but where the vector of unknowns (X1,…,X2n + 3) now contain Jx, Jy, Jz and the n top (hit, i = 1,n) and the n bottom (hib, i = 1,n) depths of the n prisms. Again, for a number mT of magnetic measurements, we have a system of mT equations with 2n + 3 unknowns. This system can be written as (9)
Considering that equations (8) and (9) have in common the top and bottom depths of the prisms, they can be rearranged as (10) where A1 is an mT ∼ 3 matrix that contains the partial derivatives of T with respect to the components of the magnetization vector, B1 is an mg ∼ 3 matrix with the derivatives of gravity with respect to the three coefficients of the contrast density function, A2 is an mT ∼ 2n matrix with the derivatives of T with respect to the top and bottom depths of the prisms and, finally, B2 is an mg ∼ 2n matrix with the derivatives of gravity with respect to the same depths of the prisms.
To solve equation (10), we follow Gill et al. (1986). We minimize the quadratic norm of residuals by means of quadratic programming, as follows: (11) subject to (12) where (Xl, Xu) are the lower and upper limits of the parameters. These constraints are very useful for introducing additional independent information such as well data, surface geology, seismic horizons, or even for fixing the upper limits of the prism to make them coincide with the surface topography. The first and second terms on the right side of equation (12) are the quadratic norms of the residuals of the magnetic and gravity data, respectively, weighted by their corresponding covariance matrix of uncertainties in the observations (CT and Cg). These matrices equalize the size of the residuals so that both sets of data are equally important. The third and fourth terms represent the quadratic norm of the derivatives with respect to the horizontal coordinates. Matrix Dt contains the derivatives with respect to the top depths and matrix Db those with respect to the bottom depths, each weighted by constants ρt and ρb. These matrices penalize the depth differences of contiguous prisms in order to obtain smooth solutions in the fashion of Constable et al. (1987), who penalize spatial variations of electrical conductivity. In the present application, we penalize spatial variations in shape as opposed to variations of physical properties.
Each of the derivative matrixes represents either the first or second derivative with respect to depth. For the first derivative, the operator produces the flattest model; for the second derivative, the result is the smoothest model (Parker, 1994). This kind of minimum roughness estimation translates in our case into structures of simple shapes, complicated only up to the point required by the data. It is easy to anticipate that the top depths will be better constrained by the data than the bottom depths, for which reason we included in equation (11) two different damping factors (ρ), one for each group of depths.
The horizontal derivatives are approximated, along the x-axis, as (13)
The explicit expression for the matrix D is then written as (Constable et al., 1987) (14)
The horizontal derivatives in the y direction are computed in the same way. Notice that the interprism distances Δ x and Δ y do not need to be the same.
Most of the applications of variable density models correspond to sedimentary basins where the density at the top is lower than at depth. But, in order to isolate the basin response from the bedrock, we require a density contrast function. As a result, the density contrast is larger at the top of the basin and lower at depth. This implies an extra difficulty in the estimation of the bottom depths, besides the fact that they are farther from the point of measurement than the tops. If the density contrast were reversed, the estimation of the bottom depths would be an easier task in spite of being farther from the gravity measurement. In general, there is interplay between the contrast density function and the top and bottom depths of the prisms. The relative stability of the different parameters depends on this interplay and on their intrinsic sensitivity to the measurements. In what follows, we present a simple sensitivity analysis for resolving the less obvious ranking between three density coefficients and the top depth of a prism.
Consider a prism with a contrast density function Δ ρ (z) = (∼0.49 + 0.068z ∼ 0.0028z2) g/cm3. We choose a prismatic square body with 1 km on a side, a variable top depth, and a total vertical length of 12 km as illustrated in Figure 3a. We can compute the derivative of gravity with respect to the top depth and with respect to the coefficients a, b, and c of the density function, and plot them as functions of the variable top depth of the prism. The derivatives are evaluated at 1 km off the side of the prism. The density function and the derivatives are plotted in Figures 3b and 3c, respectively. It can be observed that for all values of the top depth, the gravity measurement is most sensitive to the coefficient c and least sensitive to the top depth itself. The ranking of sensitivities suggests that in the inversion process the best estimated parameter will be c, then b, then a, and at the end ht, the top depth of the prism. In multiple prism inversion, we have seen that coefficients are estimated in the mentioned order, followed by the top depth and, logically, at the far end by the bottom depth.
The sensitivity to the bottom depth depends a great deal on the contrast density function. As stated earlier, the determination of this depth becomes more uncertain as the density of the sediments approaches that of the basement. Figure 4 illustrates this effect. In this case, we fixed the top of the prism and varied the bottom depth. We also moved the prism horizontally, computing the derivative of gravity with respect to the bottom depth at the point (x = 0, z = 0). The contours in Figure 4a correspond to a uniform density contrast of Δ ρ = ∼0.49 g/cm3 and those of Figure 4b to the contrast density function Δ ρ(z) = (∼0.49 + 0.068z ∼ 0.0028z2) g/cm3. It can be observed that the sensitivity of the measurement in the case of variable density is, for the longest prisms, about one-fifth of that with uniform density. This demonstrates that variable density puts extra pressure on gravity inversion methods when it comes to the determination of basin depths.
We describe in this section two applications of the algorithm. The first uses synthetic data, both gravity and magnetic, computed at two different levels for a relatively simple model of a sedimentary basin. The model has a flat upper topography that is kept fixed during inversion; we work in this case only with the bottom shape of the model. This example illustrates the different types of models that can be generated by choosing different forms for the matrix D in equation (11). The second application uses actual field data. These consist of a combination of land, airborne, and marine gravity and magnetic measurements taken over a basin that extends offshore. We use this example to illustrate the ability of the algorithm to estimate both the upper and lower boundaries of the basin, especially with multiple types of acquired data. For this, we take advantage of the marine character of the basin, whose upper topography is well known through independent bathymetric surveys. We estimate two models: one with the upper topography fixed to the known bathymetry, and one with the upper topography considered as unknown. In this way, we manage to compare the estimates of the algorithm with the actual topography of the basin, at least with its upper part.
The results of this application are summarized in Figures 5 and 6. The hypothetical model is represented in Figure 5 along with three estimations of the bottom topography of the model. The synthetic data, with noise added, is shown in Figure 6 along with the response of the model that best fitted the data.
The assumed 2D geometry in the example is realized simply by considering very elongated prisms along strike. The magnetic and gravity responses are computed for the two altitudes of ground level and 500 m above. The sediments are assumed to be nonmagnetic and the basement to possess a uniform magnetization. The density is assumed to increase with depth within the basin and is unchanging in the basement. The magnetization contrast is Δ J = ∼500 emu (I = 51°, D = 0°) and the density contrast as a function of depth is Δ ρ (z) = (∼1.0 + 0.1358z ∼ 0.00568z2) g/cm3. To the gravity and magnetic responses, we added random noise with standard deviations of 3% and 6% of the unperturbed values, respectively. We used 500 2D prisms, with the top depths fixed to the surface and the bottom depths free to accommodate themselves to the requirements of the data. The initial model has a uniform bottom depth of 500 m, a constant magnetization of ∼580 emu (I = 51°, D = 0°) and a uniform density contrast Δ ρ = ∼0.7 g/cm3. In the objective function of equation (11), we let ρt = 0, which has the effect of fixing the top depths to their initial value.
The hypothetical or reference model and three estimations by the algorithm are shown in Figure 5a. The fit to the data of the first- and second-derivative models converged in less than five iterations to the required level: 3% for gravity and 6% for magnetics. The third model, named minimum depth, is computed by replacing the matrix D of derivatives in equation (11) by the identity matrix. This has the effect of minimizing the norm of the vector of depths rather than the norm of their derivatives. The result is a shallow or minimum depth model that was too sensitive to the noise, particularly at shallow depths. Regarding the fit to the data, all we could do was to lower the misfit to around 20%. Figure 5b shows the recovered density contrast for the second derivative model. The estimated magnetization contrast was ∼441 emu (I = 60.7°, D = 1.4°). The comparisons between the responses of the second-derivative model and the synthetic noisy data at the two levels are shown in Figure 6.
The end result of the above tests is that the smoothest or second-derivative model is the most adequate of the three estimations. It was also with this option that we obtained the closest match to the density contrast function. The closeness of the flattest or first-derivative model to the smoothest one is somewhat artificial. It was achieved by using the same density function as obtained in the case of the smoothest mode. Letting the density function free, an even flatter model results than the one shown. So, even in this situation, the smoothest model is the best. It possesses the intrinsic asset of small spatial variations of the bottom geometry of the basin. This numerical experiment demonstrates that the norm of the first derivative is a more severe constraint than the norm of the second derivative. The former penalizes mainly the gradients at inflection points, where the magnitude of the derivative is the largest, so the resulting model is the flattest. In the case of the second derivative, the constraint applies largely to the peaks and valleys; the inflection points are now free and can take larger values of derivatives. In the present case, the two models fit the data equally well but, according to the arguments above and to the actual results, the most appropriate is the smoothest model. It is also the more stable when combined with a variable density function, a most important point in this case. From this experience, we decided to only use the smoothest models in actual observed data.
OBSERVED FIELD DATA EXAMPLE
Next, we used a combined set of observed field data (land, airborne, and marine gravity and magnetic) from the Ensenada Bay in the peninsula of Baja California. The bay is a relatively shallow sedimentary basin located between 31° 55′ and 31° 45′ North latitude, and 116° 50′ and 116° 35′ West longitude. The city of Ensenada is located along the coastline about 100 km south of the México-USA border. The sedimentary basin apparently opens to the west, flanked by the two islands of Todos Santos (Figure 7). The bathymetry indicates an average depth of 50 m, with a slightly deeper channel to the north of the islands and a very deep channel to the south, between the islands and the Punta Banda peninsula. The basement consists of Cretaceous andesites and granites overlaid by younger conglomerates and alluvium (see Figure 7).
Structurally, the basin is controlled by the active South and North Agua Blanca faults, in the Punta Banda peninsula (Quintanilla and Suarez, 1992) and by suspected northeast-southwest faults to the east. These faults, inferred by the formation of the basin, are covered by the conglomerates and alluvium but could be revealed by geophysical measurements. To this end, a number of local gravity surveys acquired for different purposes (Gonzalez-Serrano, 1977; Vazquez, 1980; Aguero, 1986; Cruz, 1986), were integrated and used as a regional survey. Land, marine, and airborne magnetic measurements were also available (Gonzalez et al., 2000), again acquired for a variety of purposes. The different data subsets and the particular geological setting called for a versatile and robust algorithm that could handle the interpretational problem at once.
The gravity and the magnetic data sets are presented in Figures 8a and 8b, respectively. The magnetic set, in fact composed of the subsets for different altitudes, is presented in a single plot. Free-air and Bouguer slab corrections were applied to the land gravity data. Offshore, the effect of the sea was subtracted by modeling using a density of 1.03 g/cm3 for salt water. The magnetic data were corrected by diurnal variation and by the geomagnetic field. For both gravity and magnetic data, a uniform regional was subtracted in order to isolate the basin response.
For the inversion, a grid of 30 ∼ 30 prisms was considered with a constant width of 1 km in both directions. The tops of the prisms were fixed to the topography or to the bathymetry, depending on their position. The bottom depths of the prisms were free only where there are sediments. In the areas where granitic or andesitic rocks are exposed, both the top and bottom depths were fixed to the topography using the equality constraints. Density and magnetization contrast were also freed. We assumed that granitic and andesitic rocks have the same density and susceptibility contrasts with respect to the basin sediments.
The initial model had a uniform bottom depth of 500 m and a density contrast of Δ ρ = (∼1.0 ∼ 0.05z ∼ 0.05z2) g/cm3, and a magnetization contrast of Δ J = ∼750 emu (I = 23°, D = 6°). After 15 minutes using a 300-MHz Pentium II computer and 11 iterations, we obtained the model that best fits the data. The gravity and the magnetic responses of this model are compared in Figure 8 with their corresponding data sets. It can be observed that there is a reasonably good agreement between the computed responses and the field data. The more dramatic misfit occurs for the magnetic case to the north of the area, where there is a contact between granitic and andesitic rocks, indicating that with all likelihood their magnetic properties are not as similar as expected. However, in the area of interest that corresponds to the sedimentary basin, the fit is a lot better for both gravity and magnetics.
The best fitting model has a density contrast of Δ ρ = (∼0.7 +0.0z + 0.0z2) g/cm3, indicating that there is practically no increase of density with depth, and a magnetization contrast of Δ J = ∼142 emu (I = 55°, D = 38°). The geometry of the model is presented in Figure 9, where we have removed sediments and left only the topography of the granitic-andesitic basement. The blue line in the figure is the present coastline, and the irregular black line marks the limit of the sedimentary cover. The area where both lines are farthest apart corresponds to the Maneadero valley, whose aquifer is extensively exploited for irrigation and domestic purposes. The model indicates a relatively deep basin with a barrier that rises close to sea level northeast of the islands.
A closer look at the model is shown in Figure 10, where we present the two cross-sections A-A′ and B-B′ indicated in Figure 9. The escarpment at both sides of the Punta Banda peninsula can be appreciated in Figure 10a which presents section A-A′. These escarpments are associated with the South and North Agua Blanca faults, the latter being partly responsible for the creation of the basin. The dashed line in the figure that closely follows the bottom of the model represents an alternative interpretation. This other model was obtained from an experiment that resulted from a hypothetical situation. We assumed that the bathymetry was unknown, and attempted its recovery at the same time as the bottom of the basin by letting free both the top and bottom depths of the prisms and the contrast density function. In this case, we used only the gravity data. We made this experiment to test the robustness of the algorithm when both upper and lower surfaces are unknown.
The results for section B-B′ are presented in Figure 10b. The prediction of the bathymetry from the gravity data is remarkable. Also, the three interpreted faults in the Maneadero basin are consistent in both models. These step faults complete the structure that gave origin to the basin. Notice also that the total depth of the basin of around 1500 m is the same in both models.
Our algorithm fulfills several existing needs for the interpretation of potential field data. The possibility of inverting data for different altitudes and the simplicity of including a priori information stand out as ways to ease the burden of re-leveling and of fixing parameters, respectively. On the other hand, the success of second derivatives for regularization and the stable recovery of variable density indicates that, however close, we have not reached the point of too many variables, and that the full value of potential field data is still to be realized. The work described here goes in this direction. Our algorithm handles 3D structures whose upper and lower boundaries can be arbitrary functions of position, it allows the inclusion of variable density contrasts with depth and, as input, it accepts data sets taken at different altitudes. The application to Ensenada Bay demonstrates that the algorithm can handle complex 3D models with variable density and composite data sets at different altitudes, something that falls beyond the capabilities of existing methods. We determined the geometry of the basin and identified three previously unknown step faults hidden in the sediments.
We thank CICESE and CONACYT for the financial support. Thanks also to J. M. Espinosa, R. Vázquez, A. González and CRM (Consejo de Recursos Minerales) for providing the field data sets. We also thank the anonymous reviewers and Associate Editor David A. Chapin for their valuable comments. The funding of Conacyt to project 28318-T to collect the latest marine data is also acknowledged.
Equation (3) simplifies to (A-1)
The corresponding derivatives of the vertical component of gravity with respect to a,b,c are (A-2) (A-3) (A-4)
If we integrate equation (1) over x and y, we obtain (A-5)
Deriving equation (A-5) with respect to the top and bottom depths, and using equation (2), we get (A-6)
By the fundamental theorem of calculus, we get (A-7)
A similar equation can be obtained for the ht derivative, but with the sign changed.
Departing from equation (5), we get (A-8)
The derivatives with respect the top and bottom depths can be obtained by deriving equation (5) with respect to z and then evaluating at hb or ht. That is, (A-9)
Deriving every term, we get (A-10)
A similar equation to equation (A-10) can be obtained for the derivative of the top depth. Again, the new expression has the sign changed.
- Received September 11, 2000.
- Revision received October 14, 2002.