## Abstract

The emergence of structure in reactive geofluid systems is of current interest. In geofluid systems, the fluids are supported by a porous medium whose physical and chemical properties may vary in space and time, sometimes sharply, and which may also evolve in reaction with the local fluids. Geofluids may also experience pressure and temperature conditions within the porous medium that drive their momentum relations beyond the normal Darcy regime. Furthermore, natural geofluid systems may experience forcings that are periodic in nature, or at least episodic. The combination of transient forcing, near-critical fluid dynamics and heterogeneous porous media yields a rich array of emergent geofluid phenomena that are only now beginning to be understood. One of the barriers to forward analysis in these geofluid systems is the problem of data scarcity. It is most often the case that fluid properties are reasonably well known, but that data on porous medium properties are measured with much less precision and spatial density. It is common to seek to perform an estimation of the porous medium properties by an inverse approach, that is, by expressing porous medium properties in terms of observed fluid characteristics. In this paper, we move toward such an inversion for the case of a generalized geofluid momentum equation in the context of time-periodic boundary conditions. We show that the generalized momentum equation results in frequency-domain responses that are governed by a second-order equation which is amenable to numerical solution. A stochastic perturbation approach demonstrates that frequency-domain responses of the fluids migrating in heterogeneous domains have spatial spectral densities that can be expressed in terms of the spectral densities of porous media properties.

## 1. Introduction

The migration of geofluids in the porous subsurface is of fundamental scientific importance and of considerable economic and environmental relevance to society. Many basic principles underlying flow and transport in porous media have been established, allowing widespread mineral and water resource exploitation and engineering in the regolith. However, as the world population grows and resource exploitation continues, prospective resources are becoming scarce. Increasingly, geofluid theories are being re-examined to improve the efficiency of characterization and prediction of the scale and quality of subsurface resources. One of the key recent advances is the explicit coupling of fluid dynamics with advanced geodynamic equations of state and with sophisticated biogeochemical reaction systems. This has allowed researchers to grasp the complexity of fluid–rock interactions and to note emergence of structure as a critical dynamic. Examples of emergence may include the creation of orebodies through promoted geochemical reactions, and oscillating opening and sealing of fault apertures through cyclic pumping of reactive fluids.

The complexity that gives rise to structural emergence can be embodied within the governing mathematical operators describing the physical processes, or through geometric properties of the system of interest. In seeking to establish organizational principles, it is important to understand relationships between patterns in the system. In this paper, we take a geometric view of a porous medium and seek to understand how spatial patterns of the porous medium properties will influence the spatial and temporal characteristics of a fluid filtering through under the influence of periodic boundary conditions. Our thesis is that emergence in natural systems is a function of both the nature of the local geofluid processes and of the structure of the supporting medium, and that in order to understand the generation of emergence, it will often be important to understand both the processes and the structure.

As an initial step to characterizing geometric influences on fluid migration, we consider the filtration of a fluid through a heterogeneous porous medium. To provide as much generality to the fluid process description as possible, we first identify a generalized constitutive relation for the fluid filtration and from it derive a governing equation for the propagation of boundary modes through the non-deforming porous medium. Typically, spatial data on the properties of the porous medium are scarce, which limits the ability to build comprehensive geofluid models. Spectral techniques are common means of correlating patterns at like scales, and are thus useful in inferring structure between correlated signals where data coverage is uneven. Here, a stochastic spatial structure is employed for the properties of the porous medium and a perturbation approach is used to derive a relationship between the spatial spectra of the porous medium properties and the fluid dependent variable. As shown, there is good prospect that this relationship is invertible, allowing the geometric structure of the porous medium to be estimated from observations of scalar fluid properties.

## 2. Generalized constitutive equation

We consider a non-turbulent fluid undergoing oscillating motion within a porous medium, and seek to generalize the fluid constitutive relation in order to capture a range of phenomena that may potentially apply in geofluid systems. In the simplest limit, the Darcy equation (Darcy 1856; Freeze 1994) relates the fluid flux vector, ** J**, to the local pressure gradient, ∇

*P*, 2.1that is, the flow is laminar and the fluid flux is proportional to the fluid pressure gradient. Note that, for large pressure gradients, inertial forces may outweigh viscous drag on pore boundaries and the flow may become turbulent (e.g. Furbish 1996). This is most often the case for large pore diameters (secondary porosity features and large aperture faults) where the Darcy equation no longer applies. Equation (2.1) is mathematically similar to Fick’s law for diffusing solutes and Fourier’s law for heat conduction, so, in the spirit of generalization, we will replace the pressure

*P*by the scalar dependent variable

*u*, and refer to

*κ*simply as the permeability and

*η*as viscosity. The permeability is a property of the porous medium, while the viscosity is a fluid property. Equation (2.1) can be combined with the linear continuity equation 2.2where

*λ*is a volumetric capacity coefficient, to yield the usual dynamical equation 2.3

Extensions to the linear constitutive relation (2.1) are common. Brinkman (1947) proposed the addition of a viscous drag term to model frictional forces acting on a steady fluid at grain boundaries
2.4Here, the non-negative quantity *β* is defined by *β*^{2}≡*η*′/*η*, the ratio of the effective viscosity averaged over the local flow within the pore space (*η*′) to the bulk fluid viscosity. In the absence of further information regarding grain-boundary friction averaging, *β*^{2} is often assumed to be unity, although Brinkman (1947) provides other possible estimates. As discussed by Capuani *et al*. (2003), this term is only justified for media with porosity close to unity, but has been shown to be useful for describing steady viscous flows in inhomogeneous systems, including biological tissues (Khaled & Vafai 2003).

Temporal modifications to the fluid constitutive relation are also of interest, since they allow a greater range of time-dependent phenomena to be modelled than are supported by equation (2.1). For example, the diffusive dynamical equation (2.3) admits unphysical infinite propagation of disturbances (Compte & Metzler 1997), a situation for which Cattaneo (1958) proposed the constitutive relation
2.5The non-negative quantity *τ* is a *relaxation* parameter that governs the temporal equilibration interval for disturbances to the flux ** J**. A similar relation was obtained by Maxwell (1867) in connection with the kinetic theory of gases, and was elaborated by Grad (1958). This relation, combined with the continuity equation (2.2), yields the well-known dissipating wave equation or telegrapher’s equation
2.6

Through the hyperbolic form of equation (2.6), the effect of the temporal derivative of equation (2.5) is to ensure that wave-like disturbances in the flow field propagate with velocity (King *et al.* 1998). The dissipating wave equation is widely used to describe finite correlation-time phenomena in diverse areas of physics, including granular system stress dynamics (Scott *et al.* 1998), the so-called second sound of thermo-acoustics (Chester 1963), heat conduction and thermo-elasticity (Marín *et al.* 2002; Khisaeva & Ostoja-Starzewski 2006) and in chemical reaction–diffusion (Manne *et al.* 2000).

Other temporal models are possible. Khuzhayorov *et al.* (2000) noted the phenomenological filtration model
2.7which incorporates a non-negative *retardation* parameter *ε* to modulate temporal changes in the state variable *u*. Such a constitutive relation was used by Tan (2006) to describe overshoot phenomena in viscoelastic fluids, and Wang & Tan (2008) used the same relation to perform a stability analysis for double-diffusive thermal convection of a Maxwell fluid. In analogy with the constitutive relation for an Oldroyd-B fluid, Alishaev & Mirzadjanadze (1975) introduced the dual temporal gradient relation
2.8which has since been widely applied in rheology, geophysics and in industrial applications (Tan & Masuoka 2005). Khuzhayorov *et al.* (2000) apply an upscaling technique to derive a macroscopic filtration law valid for low Reynolds and Deborah numbers that incorporates Oldroyd-B phenomena as a special case. They question the viability of equation (2.8) to model the filtration of Oldroyd-B fluids, since the predictions of their macroscopic law depart strongly from those of the phenomenological model (2.8). Nevertheless, equation (2.8) continues to be used to model viscoelastic convection phenomena in porous media, with the convective onset properties now well understood (Kim *et al.* 2003). More recent developments include the use of fractional derivative forms of equation (2.8) to model transient stress relationships in generalized Burgers’ fluids (Khan & Hayat 2008; Xue *et al.* 2008), and the potential for convective instability in doubly diffusive systems exposed to constant fluid throughflow (Shivakumara & Sureshkumar 2008).

In pursuing a simplified filtration model of viscoelastic flow in porous media, we propose the following generalized relationship for a single fluid-state variable *u* subject to Brinkman viscous drag (enhanced dissipation) and dual temporal relaxation–retardation phenomena:
2.9Here, we have absorbed viscosity *η* into the definition of *κ*. The hope is that this constitutive relation is flexible enough to describe a range of non-Darcian flow phenomena in heterogeneous porous media, but remains simple enough to yield analytical estimates of stochastic relationships in the spatial and temporal frequency domains. As we shall see, the linear nature of equation (2.9) is sufficient to ensure that such spectral relationships can be obtained in suitable limits. Chief among the limiting assumptions is the neglect of any explicit constitutive relation between the fluid stress and strain tensors, i.e. the filtration approximation. In justifying this assumption, we appeal to the limit of laminar Darcy flow, that is, where the effects of stress–strain relationships are small and can be modelled through a simple momentum law such as equation (2.9). In making this assumption, we expressly rule out turbulent effects associated with high-velocity flow through large voids. For broader application, it is understood that the general scalar state variable *u* may variously be associated with fluid pressure, concentration or temperature, according to practical context.

## 3. Modal dynamics in heterogeneous systems

Geological porous media display heterogeneity. In our model, we will seek to represent heterogeneity by allowing our material property parameters to vary in space. In particular, we deem the permeability *κ*, relaxation *τ* and retardation *ε* to take the form of random functions of location. We also assume that the Brinkman term is such that the fully heterogeneous Brinkman form may be approximated as ∇⋅(*β*^{2}*κ*(*x*)∇** J**)≈

*β*

^{2}

*κ*(

*x*)∇⋅∇

**, i.e. Brinkman drag is modelled as a homogeneous correction term. Furthermore, we assume the capacity coefficient**

*J**λ*to be uniform and homogeneous. Spatial heterogeneity of the Brinkman and capacity coefficients may be included at the expense of greater algebraic complexity. Derivation of the dynamical equation follows the usual approach of substitution of equation (2.9) in the continuity equation (2.2), yielding 3.1Since we are most interested in oscillating flows, we transform to frequency space using the Fourier transformation rule 3.2where

*h*(

*t*) is required to be integrable over time. It then follows that and . Transforming equations (2.2), (2.9) and (3.1), and taking repeated spatial derivatives of the continuity equation to eliminate

**and its derivatives from equation (3.1) eventually yields the following governing equation for responding to the frequency mode**

*J**ω*: 3.3This hyperbolic equation describes how the Fourier coefficient of

*u*for frequency

*ω*varies spatially within the porous medium in response to spatial variations in

*κ*,

*τ*and

*ε*. In the diffusive limit

*β*

^{2}=

*τ*=

*ε*=0, the governing equation reduces to the usual Helmholtz form for damped oscillations , which is well studied for oscillating modes in groundwater hydrology and heat conduction. If the relaxation and retardation fields have zero gradients, the governing equation can be written , where

*k*can be identified as the complex wavenumber of the oscillation, satisfying . The real and imaginary parts of

*k*display most curvature for small

*ω*. Since the real parts of the solutions of (3.3) are oscillatory, it is usually helpful to classify their behaviour in terms of oscillation amplitude, , and phase . Townley (1995) provides a useful mapping of the phase into units of oscillation periods, denoted Φ, 3.4The governing equation (3.3) may be integrated analytically for special classes of the spatial variations of

*κ*,

*τ*and

*ε*, and is amenable to numerical solution by standard finite-difference techniques. Appendix A discusses analytical and numerical solution of equation (3.3) in more detail. Figure 1 shows that a simple central-difference approach provides useful approximations to a wide range of dissipating wave phenomena generated by equation (3.3), and is sufficient to yield near-second-order convergence for underdamped waves with a spatially variable relaxation function

*τ*.

## 4. Spatial correlation in the weakly hyperbolic model

Having derived the equation governing the propagation of frequency modes through our generalized fluid, we are now in a position to consider how the spatial statistics of the porous medium may influence the properties of the propagating modes. One approach to this would be to simply generate a set of realizations of *κ*, *τ* and *ε* fields (according to specified internal correlation structures), run them through a numerical solver to generate the corresponding solutions and to cross correlate these against the input fields. This is a computationally expensive approach. Instead, we seek to derive, from first principles, a statistical correlation between the input and output fields that is independent of the internal correlative structures of the input fields. In order to do so, we will have to make a range of assumptions regarding the strength of the relaxation/retardation and dissipative terms in equation (3.3).

We commence by regarding *κ*, *τ*, *ε* and as spatially random fields that are stationary, that is, their means are well defined over the domain of interest with no drift terms. Consequently, we can decompose the property and solution fields into means and fluctuations (*δ*) as follows:
4.1

In equation (4.1), we have made a log transformation for all non-negative parameters. Denoting the geometric mean by the subscript ‘G’, , and , the decompositions can be expressed as 4.2

Since the property fields are now random variables, equation (3.3) requires some rationalization in order to proceed with the analysis. In particular, we note the denominators (i−*ωτ*) and (1+i*εω*+i*β*^{2}*λω*) and introduce the concept of *weakly hyperbolic* dynamics, that is where *ωτ*≪1 (weak temporal relaxation) and where both *ωε*≪1 and *β*^{2}*λω*≪1 (weak temporal retardation and weak Brinkman dissipation, respectively). In these domains, the denominators can be expanded to first order, yielding the weakly hyperbolic form of equation (3.3)
4.3Substituting equation (4.2) and neglecting gradients of the field means (a stationarity assumption) yields
4.4We now assume that the spatial fluctuations are small (a low-variance assumption), so that the exponential terms can be linearized,
4.5Substituting equation (4.5) into equation (4.4), discarding terms containing second-order and higher products of fluctuations, and noting the mean solution satisfies
4.6leads eventually to
4.7We can identify further terms that are small. By the weakly hyperbolic and low-variance assumptions, the term and can be neglected. Also, can be shown to be similarly small and can be neglected. With these reductions and defining the convenient parameter *α*=i+*β*^{2}*λω*+*ωε*_{G}, we obtain
4.8The spatial spectral characteristics of the property fluctuations must be resolved. We pursue the following Fourier transformation into reciprocal space measured by the wavevector ** k**:
4.9Equation (4.8) contains several terms that are products of space-dependent functions. In order to pursue spectral integration, we need to consider these forms explicitly. First, we consider the wavenumber transform of the term . Integration by parts yields
4.10We can reduce equation (4.10) by noting that as for all dissipative processes, so the first term in parentheses vanishes. Secondly, we are most interested in observed spectral properties near a forcing boundary, i.e. where the signal to noise ratio will be high enough to discriminate property fluctuations. With this proximity assumption, we may approximate the gradient of the mean response by a constant, i.e. . Hence, it is easy to show that
4.11where d

**Z**

_{δF}is the Fourier–Stieltjes increment (defined with support

**) of the random process**

*k**δF*(defined with support

**). Parra**

*x**et al.*(1999) discuss application of Fourier–Stieltjes forms to acoustic-wave propagation in heterogeneous systems. Applying the proximity limit equation (4.11) to the remainder of equation (4.8) and rearranging yields the following spectral relationship, expressed in Fourier–Stieltjes notation: 4.12

The Wiener–Khinchine theorem expresses the spectral density **S**_{hh} associated with the Fourier–Stieltjes increment d**Z**_{h} of a function *h*,
4.13Thus, equation (4.12) readily admits the following expression for the spectral density of :
4.14where *ς*≡*β*^{2}*λ*+*ε*_{G}.

Equation (4.14) provides a first-order expansion of the spatial spectral density of the modal fluctuation in terms of the spatial cross-spectral densities of the property fields *κ*, *τ* and *ε*. The full equation is complicated, but admits some special cases. In the case of uncorrelated property fields (which may not apply in many practical situations) the cross-spectral densities have zero means, yielding
4.15Simple diffusion/conduction phenomena, e.g. Darcy flows or heat conduction, can be recovered in the limit *ε*_{G}=*τ*_{G}=*β*^{2}=0,
4.16while Darcy–Brinkman flows can be recovered in the limit *ε*_{G}=*τ*_{G}=0,
4.17It is not until retardation and/or relaxation phenomena are enabled that cross-spectral terms appear in the modal spectral density (4.14). It is also notable that equation (4.14) tends to zero with ** k**, that is, the perturbative approximation is unable to resolve long-wavelength structures. This constraint arises from the earlier spatial stationarity assumption.

## 5. Robustness of the proximity assumption

The usefulness of the spectral density relationships derived in the previous section hinges on the robustness of the proximity assumption discussed prior to equation (4.11), i.e. the validity of the approximation . The potential usefulness of this assumption is promoted by the weakly hyperbolic approximation that renders the system more strongly dissipative than wave-like, thus reducing the occurrence of standing wave patterns that may generate strongly varying gradients in the mean response. In the following, we demonstrate that a mean value estimate of is sufficient to provide useful agreement between the input permeability spectrum and the observed modal fluctuation for systems corresponding to equation (4.16). The demonstration is performed by calculating the modal fluctuation directly from a numerical solution evaluated from a randomly heterogeneous input property field, and then comparing the spectra of the solution fluctuation and the property field, according to equation (4.16).

In practical application, the mean values of the porous medium properties will probably be determined by a routine estimation procedure using observations of the scalar dependent variable of interest. For a Darcy flow, the porous medium property is *κ*. In order to improve statistics in these small numerical experiments, we consider a dual-boundary problem over a unit square in two dimensions. By imposing synchronous Dirichlet conditions at the *x*=0 and *Λ* boundaries, and zero-gradient conditions on the other two boundaries, we can generate a system whose mean solution is given by
5.1with .

As an example of a spatially heterogeneous *κ* field, we generate a log-normally distributed spatial random field according to the algorithm of Ruan & McLaughlin (1998), and impose an isotropic Gaussian spatial autocorrelation spectrum
5.2In equation (5.2), *ρ* is the spatial integral scale and is the variance of . In order to perform the Darcian numerical solution over the unit square, we choose the following parameter values: *Λ*_{x}=*Λ*_{y}=1, *ω*=6*π*, *β*=0, *τ*=0, *ε*=0, *λ*=1 and *κ*_{G}=1/3. The spatial discretization is set at (*N*_{x},*N*_{y})=(92,92), with integral scale *ρ* set at two nodal increments *ρ*=2*Λ*_{x}/(*N*_{x}−1)≈0.022, with the corresponding wavenumber *k*_{ρ}=2*π*/*ρ*≈286. Such a small grid will not allow us to perform meaningful calculations of statistical stationarity of the property and solution fields, but it will provide a convenient test of the performance of the estimator (4.16).

To ensure conformity with the solution (5.1), we apply zero-gradient boundary conditions on both *y* edges; on the *x*=0 and *Λ*_{x} edges. Unit Dirichlet fluctuations are imposed (in phase), i.e. . Thus, the numerical solution along any line of constant *y* will approximate the one-dimensional solution (5.1). It remains to specify the variance to complete the numerical solution. Figure 2 presents (*a*) one realization of *F* with and (*b*) its spectral density. The spectral density line was averaged over *N*_{x}/2=46 randomly chosen transects through the two-dimensional density function, and is compared with the imposed Gaussian form shown by the solid line, evaluated using *ρ*=2.2*Λ*_{x}/(*N*_{x}−1) and rescaled to match at *k*=0. The 10 per cent inaccuracy in integral scale arises from the finite spectral sample offered by the discretized grid. The Nyquist wavenumber *N*_{k}=2*π*(*N*_{x}−1)/(2*Λ*_{x}) marks the maximum wavenumber generated in the spectral density and is equal to *k*_{ρ} in this simulation. The Darcian numerical solution is presented as shaded maps of magnitude and phase Φ in figure 2*c*,*d*, respectively. Dark shades indicate values near zero, while light shades indicate values near one (see scale panel). Comparison plots between the mean solution (5.1) and along the indicated line of constant *y* are shown in figure 2*e*. The choice of line was arbitrary. The final plot in figure 2*f* shows the real and imaginary parts of the fluctuation along the same line of constant *y*. At this low variance, the mean solution is a reasonably good predictor of the true response (figure 1*e*), although the fluctuation has appreciable magnitude throughout the domain (figure 1*f*). Figure 3 shows corresponding results for a realization of *F* with and integral scale *ρ*=4*Λ*_{x}/(*N*_{x}−1)≈0.044, with the corresponding wavenumber *k*_{ρ}=2*π*/*ρ*≈143, which is one half of the grid Nyquist wavenumber. The Gaussian spectrum in figure 3*b* was evaluated with *ρ*=4.4*Λ*_{x}/(*N*_{x}−1).

We are now in a position to test whether the spatial statistics of the calculated fluctuation are indeed related to those of according to equation (4.16). In order to do so, we must first estimate the mean value of . Using equation (5.1) and applying the boundary conditions , we find that the mean value of the gradient over the interval *x*∈[0,*ξ*], where *ξ*<*Λ*_{x}, is
5.3Inserting the simulated values for *Λ*_{x}, *λ*, *ω* and *κ*_{G} and trialling *ξ*=*Λ*_{x}/2 yields . Figure 4 shows the performance of equation (4.16) as an estimator of for fields, with the two different choices of variance and integral scale shown in figures 2 and 3. All spectral density lines shown were averaged over *N*_{x}/2 randomly chosen transects through the full two-dimensional density. With the above mean gradient evaluation, equation (4.16) provides a reasonable estimate of (in logarithmic terms) over a sizeable interval of wavenumbers for both integral scales. At higher wavenumbers, the estimator undershoots so that very-short-range fluctuations are not well resolved. If the log-transformation of the spectral densities is acceptable, the result is not unduly sensitive to choice of mean gradient. Variations in of ± 30 per cent yield the dashed lines in figure 4, hence there is some leeway in practical application of the proximity assumption. Increasing the variance of the input field to produces a similar quality of estimation for (not shown). This seems to be at odds with the understanding that the first-order perturbation expansion hinges on low-variance fields. However, from figure 4 it is clear that, for these single-scale Gaussian-correlated fields, the spectral estimator works best for wavenumbers less than the integral-scale wavenumber *k*_{ρ}, hence, the estimator is unable to resolve the short-range correlations (at scales less than *ρ*), where most of the high variance nature of is expressed. Another way to say this is that the governing dissipative equation upon which equation (4.16) is based is a spatial averaging operator that inherently diminishes higher wavenumber fluctuations.

## 6. Conclusion

This paper has presented governing equations for oscillating flows in porous media for the case of a generalized linear fluid constitutive equation that models viscoelastic dissipation and temporal relaxation/retardation processes. Such equations have been shown by other researchers to admit dissipating wave and overshoot phenomena, which may be of interest in the analysis of oscillating scalar dependent variables (e.g. heat, pressure, concentration) in the regolith. Simple diffusive oscillations can be recovered by appropriate choice of constitutive parameters. Turbulent processes are not considered. The linearity of the governing equations permits frequency decomposition, so that boundary forcing modes can be treated separately. A standard finite-difference numerical scheme provided good solution convergence in both diffusive and standing wave limits of the constitutive relation.

Most porous media display spatial heterogeneity in their physical, chemical and biological properties. This property richness means that forward prediction of fluid transport and reaction is problematic unless some basic spatial statistics of the property fields can be measured or inferred. Here, we have shown that, for cases where the relaxation and retardation terms are small, the spatial spectral density of oscillations in the scalar dependent variable (induced by boundary forcings) can be related to the spatial spectral densities of the porous medium property fields. This relationship was derived by a first-order stochastic perturbation approach that limits the application of the result to low-variance property fields and to measurements taken close to the forcing boundary. The usefulness of the spectral relationship was examined for a Darcy flow problem over a randomly heterogeneous permeability field. The spectral relationship was found to be a good estimator at low and moderate permeability variances at least, and was not unduly sensitive to the spatial average of mean gradients near the forcing boundaries (a required parameter). Exhaustive testing of the spectral relationship for heterogeneous relaxation and retardation property fields is beyond the scope of this paper, but is relevant to the study of strongly dissipating wave-like phenomena in granular and sedimentary systems.

The value of the spectral relationship derived here is that it is invertible for a useful interval of wavenumbers. That is, it may be used to infer spectral characteristics of the porous medium property fields over a range of spatial scales from observations of the space–time spectra of the scalar dependent variable. The predictive power of the estimator is least in the small- and large-wavenumber limits, and greatest in wavenumber ranges corresponding to distances of several correlation scales of the input property fields. Typically, *in situ* observations of pressure, heat or concentration are inexpensive in comparison to *in situ* measurements of porous medium properties, so the inversion offers potential benefits in terms of property characterization. However, in order to provide robust estimates of the spatial correlation structure of the problem, it seems to be necessary to extend the accuracy of the estimator to higher wavenumbers. Furthermore, optimal sampling techniques may need to be considered in order to reduce uncertainty of property estimates based on small numbers of observations of the scalar dependent variable, which is most likely to be the case in practical application.

## Appendix A. Analytical and numerical solutions

In the following text, we compile some useful analytical solutions to the governing equation (3.3). Consider the diffusive limit *β*^{2}=*τ*=*ε*=0 over the finite domain *x*∈[0,*Λ*], with a Dirichlet condition at the origin, i.e. , and a Neumann condition at the upper boundary, i.e. . The simple case where *κ* is a constant yields
A 1with . The permittivity of the medium is parameterized by the modulus of *φΛ*: where |*φΛ*| is large (e.g. high frequencies), the migrating disturbance is rapidly dissipated, and where |*φΛ*| is small (e.g. low frequencies), the migrating disturbance can persist over long distances. For the case *κ*=*κ*_{0}+*κ*_{1}*x* and Dirichlet conditions at both ends of the domain and , equation (3.3) is also soluble, yielding
A 2where *K*_{0} is the modified Bessel function of the second kind of order zero, is the regularized confluent hypergeometric function
A 3and the auxiliary quantities *φ*_{0} and *φ*_{1} are given by and , respectively.

Now consider a system where constant values of *κ* and *τ* apply throughout the domain, i.e. a constant value of temporal relaxation is operative, and where Dirichlet boundaries apply at each end, as before. Retardation phenomena are neglected. The following solution is obtained for non-zero Brinkman dissipation *β*:
A 4where the scaling constant *γ* is given by
A 5This solution represents a dissipating mode penetrating simultaneously from each boundary of a uniform medium. Relative phase differences between the boundary conditions (i.e. asynchronous forcings) can be adjusted by specifying and in complex form. A complementary case is presented by a medium with constant values of *κ* and *ε*, but where relaxation phenomena are neglected. In solving equation (3.3), we obtain (A4) again but with a different definition of the scaling constant,
A 6Finally, consider a system where constant values of *κ* and *λ* apply throughout the domain, where Dirichlet boundaries apply at each end, and where *τ*=*τ*_{0}+*τ*_{1}*x*, i.e. the relaxation time, is a linear function of position. In the absence of Brinkman dissipation, the following solution is obtained:
A 7where the following auxiliary identifications are made: , *σ*_{0}=(−i+*τ*_{0}*ω*), *σ*_{1}=(−i+*τ*_{0}*ω*+*Λτ*_{1}*ω*). Solutions (A1), (A2) and (A5)–(A7) are presented in figure 1*a*–*e*, respectively. All results shown in *a*–*e* of figure 1 were calculated with (*N*_{x},*N*_{y})=(51,5) and with parameter values expressed in arbitrary units. Figure 1*a*,*b* shows diffusive solutions for equations (A1) and (A2), respectively, with common parameter values *ω*=8*π*, *λ*=1 and *β*=*τ*=*ε*=0. Figure 1*a* uses a unit Dirichlet boundary at *x*=0, a zero-gradient boundary at *x*=*Λ* and constant *κ*=0.2. Figure 1*b* uses unit Dirichlet boundaries at *x*=0,*Λ*_{x} (also applied in figure (*c*–*e*)), *κ*_{0}=0.2 and *κ*_{1}=0.5. Figure 1*c*,*d* shows hyperbolic solutions for (A4), (A5) and (A6), respectively, with *ω*=4*π*. Figure 1*c* uses constant *κ*=0.2, *λ*=1, *τ*=0.15 and *β*=*ε*=0, and figure 1*d* uses constant *κ*=0.2, *λ*=1, *ε*=0.3 and *β*=*τ*=0. Figure 1*e* shows solutions for the variable hyperbolic case (A7) with *τ*_{0}=0.1,*τ*_{1}=0.2,*ω*=4*π* and constant parameters *κ*=0.35, *λ*=1 and *β*=*ε*=0.

Direct integration of the generalized governing equation can be achieved via finite differences as follows. Consider a two-dimensional domain (*x*,*y*)∈[0,*Λ*_{x}]×[0,*Λ*_{y}], divided into a regular rectangular grid of *N*_{x}×*N*_{y} nodes, over which the numerical solution *u*_{i,j} is defined (extensions to three dimensions are straightforward). Coordinate mappings are given by (*x*_{i},*y*_{j})=((*i*−1)Δ*x*,(*j*−1)Δ*y*) for *i*=1,2,…,*N*_{x} and *j*=1,2,…,*N*_{y}, with Δ*x*=*Λ*_{x}/(*N*_{x}−1) and Δ*y*=*Λ*_{y}/(*N*_{y}−1). If we assume that *κ* is defined and piecewise continuous at all points of the domain, then we arrive at the following derivative discretizations:
A 8and
A 9In (A8) and (A9), the right-hand sides are expressed in terms of relative coordinate indices, so that , etc. Assembling all the terms in equation (3.3) yields
A 10If we view only the nodal values *κ*_{i,j} as being known, then using the normal effective diffusion result for slabs of differing *κ* in series (Crank 1975), mid-nodal values may be estimated using harmonic means,
A 11

As shown in figure 1, the stencil defined by (A9) and (A10) provides useful accuracy for simple and mildly heterogeneous media. For very highly heterogeneous media, more sophisticated discretizations are possible (Alcouffe *et al.* 1981; Shashkov & Steinberg 1996), but are not discussed further here.

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

- © 2010 The Royal Society