Fluid Dynamics Research 40 (2008) 273 – 309
Dispersion modeling by kinematic simulation: Cloud dispersion model J.C.H. Funga,∗ , R.J. Perkinsb a Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong b Laboratoire de Mécanique des Fluides et d’Acoustique, Ecole Centrale de Lyon, France
Received 22 August 2004; received in revised form 22 May 2007; accepted 22 June 2007 Available online 20 November 2007 Communicated by M. Oberlack
Abstract A new technique has been developed to compute mean and fluctuating concentrations in complex turbulent flows (tidal current near a coast and deep ocean). An initial distribution of material is discretized into any small clouds which are advected by a combination of the mean flow and large scale turbulence. The turbulence can be simulated either by kinematic simulation (KS) or direct numerical simulation. The clouds also diffuse relative to their centroids; the statistics for this are obtained from a separate calculation of the growth of individual clouds in small scale turbulence, generated by KS. The ensemble of discrete clouds is periodically re-discretized, to limit the size of the small clouds and prevent overlapping. The model is illustrated with simulations of dispersion in uniform flow, and the results are compared with analytic, steady state solutions. The aim of this study is to understand how pollutants disperses in a turbulent flow through a numerical simulation of fluid particle motion in a random flow field generated by Fourier modes. Although this homogeneous turbulent is rather a “simple” flow, it represents a building block toward understanding pollutant dispersion in more complex flow. The results presented here are preliminary in nature, but we expect that similar qualitative results should be observed in a genuine turbulent flow. © 2007 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. PACS: 47.27.Eq; 92.10.Lq; 92.20.Ny Keywords: Turbulence dispersion; Kinematic simulation; Concentration fluctuation; Cloud dispersion model
∗ Corresponding author. Fax: +852 2358 1643.
E-mail address:
[email protected] (J.C.H. Fung). 0169-5983/$32.00 © 2007 The Japan Society of Fluid Mechanics and Elsevier B.V. All rights reserved. doi:10.1016/j.fluiddyn.2007.06.005
274
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
1. Introduction There are many environmental problems involving the release and subsequent dispersion of material, either in the atmosphere or in rivers and oceans. Often it is necessary to be able to compute the concentrations which will occur, and a number of techniques have been developed to do this. Most of these compute only the time-mean or ensemble-mean concentrations but in many situations this is insufficient because there are processes which occur on time scales that are short compared with the time that it takes for the plume, or cloud, of matter to be dispersed. Such physico-chemical reactions or biological reactions would be quite incorrectly estimated if the lower level of mean concentration was used rather than the high concentration that occurs over short time and length scales. In these situations it is necessary to consider how material disperses in a turbulent flow. When material is released into a turbulent flow in the form of a cloud, it begins to disperse. Eddies of the size of the cloud and smaller act to deform the cloud, stretching it in some directions, compressing it in others, and drawing out ‘tendrils’ of material. Eddies of the size of the cloud and larger move the cloud around, but as the cloud spreads outward more and more of the scales of motion are involved in deforming the cloud, so deformation occurs over progressively larger length scales. However, as the cloud is drawn out and distorted the length scales become smaller. This leads to mixing with surrounding fluid to create high local concentration gradients, which enhances diffusive flux at a molecular level. These aspects can be modeled statistically (if the parameters are known), but such models do not relate mixing to the actual eddy structure. Recent research shows that structural diffusion (Fung et al., 1991) (diffusion depending on the spatial structure of the flow field) or stirring (Garrett, 1983) is also important because it leads to finer scale variations in the concentration at any instant, and hence to larger concentration fluctuations. In any single release of material the cloud does not spread evenly; some eddies draw material out in one direction, others draw material in different directions. Averaged over many releases in homogeneous, isotropic turbulence, the concentration distribution would be smooth and radially symmetric, but in any single realization it is not so (Fig. 1a). It is clear then, that the instantaneous structure of the velocity field plays an important role in determining the fluctuating concentration. Unfortunately, there are very few techniques which are capable of computing the concentration distribution in a single realization of a turbulent velocity field. Theoretical treatments of the diffusion of material in a turbulent flow have generally followed one of the four approaches—gradient transport models, stochastic models, full numerical simulations and kinematic simulations (KS). In gradient-transport models, the conservation equation for material released into an incompressible velocity field u(x) is jc/jt + u · ∇c = 0, where c is the local concentration of material (per unit volume). Writing u = (u, v, w) and c as the sum of mean and fluctuating components (e.g. c = c + c ), and then averaging jc ju c jv c jw c + u · ∇c = − + + . (1) jt jx jy jz In the gradient-transport model it is assumed that the local turbulence causes a net transfer of material down a concentration gradient, so that terms such as u c may be re-written as −Kx jc/ ¯ jx, where Kx is the diffusion coefficient in the x-direction. (The physical basis for this argument is the same as for the turbulent transport of momentum: an eddy with length scale l and characteristic velocity u transports a parcel of fluid with concentration c from x to x + l, where the concentration is c + (jc/jx)l. From this,
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
275
single realisation
ensemble mean
Fig. 1. (a) The difference between the dispersion of a cloud in a single realization of a turbulent velocity field (continuous line) and the ensemble-average of many realizations (dashed line). (b) Schematic showing all of the small scale flows, large scale flows and the mean in the simulation of KS2.
u c = −u l(jc/jx), from which u l ∼ K and u c ∼ −K jc/ ¯ jx.) Substituting this expression into (1) we obtain jc j jc + u · ∇c = Kij , (2) jt jxi jxj where Kij is the diffusion coefficient tensor to account for different transport rates in different directions (Batchelor, 1949). This expression (2) is known as the advection–diffusion equation. Analytic solutions exist for certain special conditions (e.g. u = (u, 0, 0)) but in general the advection–diffusion equation has to be solved numerically using a finite difference scheme. Typically, values for the diffusion coefficient K are found to vary with position, time of travel and the scale of the diffusion process. In stochastic models the dispersion of material is modeled by computing the trajectories of many thousands of particles passing through the source. If the material is a scalar then the particles are essentially fluid elements. The basic assumption is that the velocity of a fluid element in homogeneous isotropic turbulence can be modeled by the Langevin equation, from which we obtain t 2t u(t + t) = 1 − u(t) + u , (3) TL TL where TL is the Lagrangian integral timescale, 2u is the variance of the fluid velocities and is a vector of random numbers drawn independently from a Gaussian distribution, with zero mean and unit variance. The correlation term (1 − t/TL ) models the fact that the fluid velocity experienced by a fluid element is
276
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
not entirely random, but shows some correlation with its previous value. Particle trajectories are calculated by integrating (3) (together with any mean velocity field) for as long as necessary. Concentrations can be obtained by computing many thousands of such trajectories, all starting at the source. Dispersion in strongly inhomogeneous and skew turbulence has been modeled by Thomson (1984, 1986) and van Dop et al. (1985) using non-Gaussian distributions for the random forcing. Thomson (1987) shows that non-Gaussian distributions can lead to models which do not satisfy physical requirements such as the well-mixed condition. To overcome this, Luhar and Britter (1989) developed a model for inhomogeneous skew turbulence (a convective boundary layer) using the sum of two Gaussian distributions to generate the random forcing. The stochastic models described so far all compute the trajectories of single, independent particles. This means that they can be used to estimate mean concentrations c, but not concentration fluctuations, because the fluctuations depend on the structure of the velocity field, and this structure provides correlations between the trajectories of neighboring particles. To overcome this restriction, two-particle stochastic models have been developed (e.g. Durbin, 1980), in which pairs of particles are released from the same point and their trajectories calculated. This requires information about the spatial correlation between velocities, particularly at small separations. If at an initial time a group of fluid particles is identified by some passive scalar property (such as a weak source of heat or dye), then the disorderly turbulent motion of the fluid particles provides a direct mechanism for the subsequent mixing of the diffusing scalar property with the surroundings. In particular, the mean concentration can be expressed in terms of “single-particle” statistics in which each particle is treated as an independent member of an ensemble. Furthermore, the concentration variance is determined by “two-particle” statistics concerned with a pair of fluid particles diffusing relative to each other. The relative motion of fluid elements also governs the change in shape of an admixture cloud released in turbulent flow. The above considerations illustrate the importance of understanding the relative motions of pairs and groups of fluid particles. In order to model processes such as mixing and combustion we also need statistical information about the relative displacement of particle pairs in a turbulent flow (Batchelor, 1952). In many situations we 2 need to predict the local fluctuations c (t) not just the mean since peak concentrations can be more important. For example, if poisonous gas is released into the air, especially in a confined volume, what matters to the potential victims of such an accidental discharge is not the mean concentration of some distance from the source but the likely maximum concentration they may encounter, and clearly the 2 fluctuation c (t) is a better measure of this. Concentration variance is intimately related to the statistics 2 of the trajectory of pairs of particles, and estimates of c (t) can only be obtained from the statistics of two-particle dispersion. The need to predict values of this quantity has led to significant advances in modeling concentration fluctuations in recent years based upon two-particle stochastic models (Durbin, 1980; Sawford, 1985; Thomson, 1990). An alternative approach is to simulate the space-time evolution of individual turbulent velocity fields and then compute concentration distributions by releasing many thousands of particles into the velocity field and computing their trajectories. Advances in computing power have now made it possible to carry out direct numerical integration of the Navier–Stokes equations (DNS), but these simulations are currently limited to turbulent Reynolds numbers of about 460 or less (Vincent and Meneguzzi, 1991; Gotoh, 2002). Even then, the computations are restricted to simple boundary conditions, and require the most powerful computers and large amounts of computing time.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
277
Because of these limitations the velocity field has to be approximated. One approach is to compute solutions to the dynamical equations for large scale motions, and use a sub-grid scale model to include the effects of the small scales on the large scales (e.g. as an eddy viscosity). This is the technique used for synoptic weather forecasting, large-scale ocean modeling and large-scale turbulence modeling (large eddy simulation, LES). As with DNS, LES require large mainframe computers and large amounts of computing time (Flohr and Vassilicos, 2000). An alternative approach, based on an idea of Kraichnan (1970), is to simulate the kinematics of the velocity field without solving the dynamical equations. This technique is known as “Kinematic simulation”. The velocity field is represented as the sum of a number of random Fourier modes, in time and space. The amplitudes of the modes are chosen so that each realization of the velocity field satisfies continuity, and an ensemble-average of many realizations satisfies a prescribed energy spectrum. Typically 50–100 modes are sufficient to cover the full range of eddies in a fully turbulent flow when the Reynolds number is of order 104 (Fung et al., 1992). Although KS does not contain any dynamical effects (there is no interaction between the modes, no vortex stretching, etc.) it does simulate some of the effects of dynamical interactions, by specifying a spatial spectrum and the form of the spatial–temporal relation in each mode. Velocity fields calculated from KS have been compared with theoretical and experimental results and with velocity fields obtained from DNS and LES. The results which are particularly important for dispersion studies concern Lagrangian statistics and the structure of the velocity field. Computed Lagrangian statistics (e.g. velocity correlations, spectra of pressure fluctuations) agree well with other computed and measured data (Fung et al., 1992, 2003). The structure of the velocity field has been compared with that of velocity fields computed by DNS, by assessing the proportions of the flow occupied by different types of flow structures (swirling eddy regions, straining regions, streaming regions, etc.). In general, the agreement is good, except that the lengths of the vortical regions tend to be too small (Fung et al., 1991). When a volume of matter is released into a turbulent flow the elements within the volume move in slightly different directions so that the distance between neighboring particles on average increases. This process can be modeled using KS to generate a velocity field and then computing the trajectories of a large number of elements within that initial volume. This method has been used to compute the evolution of material surfaces in turbulent flows (Malik, 1990; Drummond and Münch, 1991). However, it is found that even when the velocity field u(x, t) is known, approximations are necessary in estimating the concentration field C(x, t) from computations of individual trajectories because the elements rapidly disperse. A number of algorithms have been proposed for re-discretizing the concentration field at different stages, to control the number of elements that have to be tracked, whilst maintaining the resolution of the concentration field. Malik et al. (1991) describe a method by which clouds of material are released and assumed to spread relative to their centroids. The velocity of the centroid of each cloud was modeled as a stochastic process, independent of all other clouds. The statistics for the growth of the clouds and the velocities of their centroids were derived from simulations of the motion of single clouds in a velocity field generated by KS. Periodically the ensemble of clouds was re-discretized, both to limit the size of the clouds and to prevent individual clouds from approaching each other too closely. The results of that model compared satisfactorily with field data and with results from another, more conventional model. The method described in this study is an adaptation of that algorithm. The continuous release of matter is still modeled by the release of discrete clouds, but the velocity of each cloud centroid is obtained directly from a KS of the velocity field. The clouds grow according to results from a small-scale simulation of
278
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
the motion of single clouds. The ensemble of discrete clouds is periodically re-discretized to prevent individual clouds from growing too large, and to maintain the same spatial resolution throughout the cloud. In Section 2 the method of the cloud dispersion model (CDM) is introduced. In Section 3, we discuss how to choose the appropriate parameters for the CDM and in Section 4, we show the result of simulations. Section 5 is the conclusion.
2. The cloud dispersion model In this section we describe the CDM and its use to simulate dispersion from a continuous release into a turbulent velocity field. We assume that the release takes place from an extended source (having an area of order L2R ), that there is a mean flow (U) and that concentrations are to be calculated on a grid with a spacing of LG . The CDM can be divided into two simulations—a small-scale simulation (KS1) and a large-scale simulation (KS2). In KS1 we simulate the growth of a single cloud from initial radius LR to a radius of order LG , using KS to generate an isotropic, homogeneous turbulent velocity field with an appropriate energy spectrum (see Section 3), and wavenumber range kG < k < kN (see Figs. 1b and 2). By taking the ensemble-average of many such simulations we compute two statistical processes—the rate at which the cloud spread increases with time, and the variation of the centroid velocity as a function of time. These statistics are then used in KS2. In the large-scale simulation the velocity field is modeled using KS, with wave numbers in the range k1 < k < kG —in other words with length scales larger than the scale of the grid on which the concentrations are to be computed. The continuous release of material is discretized into the release of “clouds” at fixed time intervals. Each cloud is represented by three parameters—the quantity of material in the cloud, the
Fig. 2. The energy spectrum showing the wavenumber ranges for KS1 and KS2: (a) for S1, (b) for S2.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
279
location of its centroid and its radius. Once a cloud has been released it is advected by the velocity field; this has three components—the mean velocity (U), a turbulent velocity field (u(x, t)) generated by KS and a random component (ur ) modeling the effects of wave numbers above kG , which are not present in the large-scale simulation. (The statistics for the random component are obtained from the results of KS1.) The cloud trajectory is obtained by a Lagrangian integration (i.e. following the cloud centroid) of the local velocity. As the cloud moves, it is also growing; the rate at which it grows is obtained from KS1. Eventually the cloud spreads to such a size that the radius is of the same size as the smallest length scales in the flow (i.e. rc (t) ∼ LG ); at this point it is no longer realistic to regard the cloud as a single entity, because the turbulent velocities are no longer correlated over the whole extent of the cloud. Thus, the cloud is split into several smaller, independent clouds. These clouds then follow the same process of advection, growth and splitting. If each cloud is split into N smaller clouds every Ts time steps then the total number of clouds at time t will be of the order of N t/Ts . For any realistic simulation the number of clouds becomes very large, very rapidly, and this causes two problems. Firstly, the number of clouds grows more rapidly than the area of the cloud (which increases as t 2 ) so that the clouds are being packed closer and closer together. Secondly, the number of clouds rapidly exceeds the practical limits imposed by computer memory. So instead of simply splitting the clouds the whole distribution of clouds is re-discretized; this involves the splitting of some clouds but the merging of others, so that the overall spatial density remains constant and the number of clouds only increases slowly (and controllably). In the following sections we describe KS, and each of the processes involved in the small-scale simulation (KS1) and the large-scale simulation (KS2).
2.1. Kinematic simulation Most simulations of turbulent flow fields use exact or approximate solutions to the equations of motion which are continuously solved to represent the evolving flow field. In a direct numerical simulation of turbulence the flow is calculated from first principles without any closure assumption. All that is needed are proper initial and boundary conditions, an accurate and efficient numerical solution scheme and a large computer (Rogallo and Moin, 1984). Although some progress has been made in the efficiency and accuracy of computational algorithms, particularly in the adaptation of spectral methods, the primary limitation on our ability to simulate high-Reynolds-number turbulence is the speed and memory size of the computing hardware. Though a turbulent Reynolds number of 150 has been achieved (Vincent and Meneguzzi, 1991), this is still too low for the study of many turbulence phenomena found in geophysical flows. This restriction on Reynolds number is avoided by computing solutions to the dynamical equations only for the large scales, and modeling the small scales in terms of a local eddy viscosity acting on the large scales. This is the method of LES. In this paper random flow fields are generated from certain statistical distributions that are known from measurements or direct simulations, but the flow fields do not necessarily satisfy the dynamical equations nor do they have all the known statistical distribution. KS is a technique for synthesizing a time-varying velocity field with a given energy spectrum. The method was originally proposed by Kraichnan (1970) and has been considerably developed since then (Turfus and Hunt, 1987; Fung et al., 1992). A full description of its use to simulate three-dimensional (3D) isotropic homogeneous turbulence is given in Fung et al. (1992), and we will only summarize two-dimensional (2D) simulations of KS. The KS of homogeneous
280
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
turbulence presented here is not grid-based, and does not require the solution of any set of equations, so it can be easily programmed on any computer, and is ideally suited to parallel computations. Most of the development and testing of KS has involved 3D velocity fields. However, the CDMs described in this study requires a 2D velocity field, and there are some important differences between 2D and 3D turbulence. The circulation of the ocean is severely constrained by density stratification. A water parcel cannot move from one surface of constant potential density to another without changing its salinity or its potential temperature. There are virtually no sources of heat outside the sunlit zone and away from the bottom where heat diffuses from the lithosphere, except for some interesting hydrothermal vents in special regions. Because of the stratification of the ocean, we will focus our study of dispersion in stratified water, where motions are constrained along isopycnal surfaces, and the flow is 2D. Hence this paper focuses on 2D flow. In this section, therefore, we provide a detailed description of KS and its use to generate 2D velocity fields (see also Vassilicos and Fung, 1995; Fung and Vassilicos, 1998). The velocity field is generated by the sum of random Fourier modes Nk u(x, t) = [An cos(n .x + n t) + Bn sin(n .x + n t)],
(4)
n=1
where Nk is the number of modes in the simulations and the Cartesian coordinates of An , Bn and kn are given by An = An (cos n , − sin n ), Bn = Bn (− cos n , sin n ) and kn = kn (sin n , cos n ). The angles n are random and uncorrelated with each other and the velocity field (4) is incompressible because An · kn = Bn · kn = 0 for all n. The random amplitudes An and Bn are chosen independently from a Gaussian distribution with zero mean and a variance (Fung et al., 1992) given by A2n = Bn2 = E(kn )kn ,
(5)
where E(k) is a prescribed Eulerian energy spectrum in the range 2/L = k1 k kNk = 2/ and such that E(k) = 0 outside this range. In this way, we obtain a particular realization of a velocity field whose spectrum is E(k). kn =(kn+1 −kn−1 )/2 for 2 Nk Nk −1, k1 =(k2 −k1 )/2 and kNk =(kNk −kNk −1 )/2. The distribution of wavenumbers kn is geometric, i.e. kn = k1 n−1 ,
(6)
where is a dimensionless number which is a function of L/ and Nk because kNk = 2/. (Hence = (L/)1/(Nk −1) .) The frequencies n in (4) determine the unsteadiness associated with wavemode n. We chose a model (Vassilicos and Fung, 1995; Fung and Vassilicos, 1998) where the unsteadiness frequency n is proportional to the eddy-turnover time of the wavemode n, i.e. n = kn3 E(kn ),
(7)
where is a dimensionless constant. The KS velocity fields simulated here are stationary in time and the auto-correlation of Lagrangian velocities following fluid elements is R L ( ) = exp(− /TL ) (see Fung and Vassilicos, 1998; Flohr and Vassilicos, 2000). However, the appearance of coherent structures can introduce long decorrelation times in the flow and produce anomalous diffusion. This behavior cannot be captured by “kinematic simulations” with finite decorrelation times.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
281
It may be worth noting that in KS we prescribe a spatio-temporal structure to the flow via incompressibility and Eqs. (4)–(7). This spatio-temporal structure incorporates eddying and straining flow regions (Fung et al., 1992; Vassilicos and Fung, 1995; Fung and Vassilicos, 1998). KS has the same advantage as spectral methods over finite-difference or finite-element methods in representing the velocity field. But in a dynamical simulation it is necessary to include modes for every wavenumber in the representation, so the number of points Ng ∝ k3 ∼ Re9/4 , where k is the highest wavenumber in the system. In KS fewer modes are needed to represent all the major kinematical features of the flow field, and it is found empirically that Nk ∝ k ∼ Re3/4 . KS has also been used to study the settling velocity of particles depending on the size of the particle and on the drag and body forces acting on it (Fung, 1993, 1998; Maxey, 1987). The transport of scalar fields by turbulent flows is an important process in many physical situations ranging from the dynamics of the atmosphere and the ocean to chemical engineering. Specific examples of scalars are provided by pollutant density, temperature or humidity fields and the concentration of chemical and biological species (Shraiman and Siggia, 2000). Issues of transport and mixing in turbulence are directly related to the properties of fluid trajectories. The problem is thus often addressed using Lagrangian techniques (Yeung, 2002). The KS provides a Lagrangian model of turbulent diffusion, based on a simplified incompressible velocity field, with a proper energy spectrum. This model reproduces very well the Lagrangian properties observed in laboratory experiments (Nicolleau and Vassilicos, 2003), as well as in DNS (Malik and Vassilicos, 1999). The computational simplicity of the KS allows us to consider very large inertial ranges: a ratio of scales of ∼ 104 is easily accessible with moderate computer resources. The KS thus turns out to be a powerful tool to study issues of dispersion in turbulent flows. 2.2. Small-scale simulation The small-scale simulation is used to compute the motion and growth of a single cloud, from its moment of release to the time when it reaches a size comparable with that of a cell in the concentration grid. We denote the length scale of the cloud by LR , and the length scale of the concentration grid by LG , and we assume that only wavenumbers in the range kG (=2/LG ) to kN (=2/LR = kR ) are involved in spreading the cloud. Eddy motions at wavenumbers below kG can advect the whole cloud, while wavenumbers above kN mix the cloud. Of course this is a simplification, because the clouds do not spread out circularly, but tend to get drawn into long thin streams with occasional bulges. However, increasing the number of wavemodes increases the computational time, and we show in Section 4.3 that truncating the wavenumber range does not affect the cloud statistics significantly over the range of cloud sizes modeled in this simulation. To check the effect of truncating the wavenumber range larger than kNk , we have also calculated 2 with the unresolved motions of turbulence for k > kNk as a stochastic forcing (as in Section 2.3.4) with a variance equal to the total velocity variance of unresolved motions. The results for the full and truncated wavenumber ranges are less than 2%, hence it is clear that the Fourier modes with scales smaller than the cloud size at release can be omitted in the simulation. We have generated many different realizations of a 2D velocity field using KS with the appropriate wavenumber range. In each realization of the velocity field a single cloud is released into the velocity field; as it is advected and deformed by the velocity field we compute the centroid velocity and the growth of the cloud as functions of time. From the ensemble-average of the results from all the smallscale simulations we obtain analytical expressions for the cloud radius and the variance of the centroid velocities, as functions of time.
282
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
The initial cloud is discretized into a large number of “particles” covering the same area as the extended source; typically in these simulations we have used square clouds composed of 40 × 40 particles. (The exact shape of the initial cloud is unimportant.) The particles are all represented by their locations xi (t). The velocity at the location of each particle is given by Eq. (4), and the particle displacements at each time step are computed using a second-order scheme due to Drummond et al. (1984). They show that this scheme has the same accuracy as a fourth order Runge–Kutta scheme, but requires fewer computations. Each particle is then advected according to (4). At selected instants the velocity of the cloud centroid and the spread of the cloud relative to its centroid are calculated using the following relationships. (a) Velocity of cloud centroid: The cloud centroid at time t is given by Np 1 xc (t) = xi (t) Np
(8)
i=1
from which the velocity of the cloud centroid is given by Np d 1 ui (t). uc (t) = xc (t) = dt Np
(9)
i=1
The large-scale simulation actually requires the variance of the centroid velocities, obtained from an ensemble-average over all the realizations. If we denote the realization number by j (j = 1 to NR ) then the mean of the centroid velocities is given by NR 1 u(t) = uc,j (t), NR
(10)
j =1
and the variance of the centroid velocities is given by 2u (t) =
NR 1 uc,j (t).uc,j (t) − u(t) · u(t). NR
(11)
j =1
(b) Cloud spreading: The radial spread of the particles relative to the cloud centroid is given by 2r (t) =
Np Np 1 1 2 |xi (t) − xc (t)| = |xi (t)|2 − |xc (t)2 |. Np Np i=1
(12)
i=1
Note that this is not the same as the spread relative to the release point (= N1p
Np
2 i=1 |xi (t)| ).
2.3. Large-scale simulation The large-scale simulation models the advection and growth of a cloud or plume of material in a turbulent field, and computes the time-dependent concentration at the nodes of a grid. For an instantaneous release, the cloud is modeled as a single cloud with an initial area equal to that of the source. For a continuous release the discharge is represented by the release of discrete clouds at fixed time intervals. Each of these clouds then has the same initial area as that of the source.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
283
Once released, the clouds grow and are advected by the velocity field. The plume is assumed to be well mixed at scales smaller than the resolution of the concentration grid, so the turbulent velocity field only contains motions at length scales greater than or equal to the grid spacing. The growth of the cloud, from its release size to the size of the grid cells, is assumed to follow the ensemble-averaged results from the small-scale simulation. The set of clouds that have grown larger than a grid cell is periodically re-discretized to create a set of smaller clouds; this involves some splitting and some merging of clouds but it generally results in an increase in the number of clouds in the simulation. 2.3.1. The large-scale velocity field As with the small-scale simulation, the velocity field is generated using a 2D KS. This time the lowest wavenumber is determined by the overall length of the simulation (LS ) and the highest wavenumber is determined by the resolution of the concentration grid (LG ), so k1 = 2/LS k 2/LG = kG (see also Fig. 1b). The energy at wavenumbers above kG (which is not included in this simulation) will increase the rate at which the plume spreads, in two ways—it will cause the individual “clouds” in the simulation to grow and it will increase the motion of those clouds. The growth of the small clouds is modeled using statistics from the small-scale simulations on the growth of clouds, and the increased motion of the clouds is incorporated into the large-scale simulation by giving each cloud centroid an extra “random” velocity at each time step. The statistics for the stochastic element in the simulation are obtained from the small-scale results for the variance of the cloud centroid velocities. The velocity field may also contain a mean flow, U. The turbulent velocity field generated by KS is computed in a frame of reference which moves with the mean flow (i.e. the turbulence is being advected by the mean flow) so the whole simulation is carried out in a frame of reference moving with the mean flow. This means that if the source is located at xR (0) at t = 0 then its location in this moving frame of reference at any subsequent time t is given by xR (t) = xR (0) − Ut.
(13)
When the concentrations are computed the simulation is returned to a stationary frame of reference by the addition of Ut to all coordinates. 2.3.2. The cloud concentration distribution Each of the discrete clouds in the simulation is represented by three parameters—the quantity of material in the cloud (Q), the position of the cloud centroid (xc , yc ) and the radius of the cloud (). The clouds are assumed to be radially symmetric about their centroids, with a Gaussian concentration distribution Q (x − xc )2 + (y − yc )2 C= , (14) exp − 22 2 2 where Q is the total amount of material in the cloud and is a measure of the spread of the cloud. Although the concentration is defined everywhere in the xy-plane, and is everywhere non-zero, the total amount of material outside a circle of radius r / approaches zero rapidly as r / increases. Making the usual substitution x = xc + r cos , y = yc + r sin and changing the variable of integration,
284
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Table 1 Ratio of material in a Gaussian concentration distribution as a function of distance from the center r Q Q
1.0 0.393
2.0 0.865
2.5 0.956
we obtain from (14): 2
2 r / 1 r 1 r 2 Q 2 = drd , r exp −
= 1 − exp − Q 2 2 22 0 0
3.0 0.989
(15)
where Q is the amount of material within a radius r /. Evaluating this for different r / we obtain the results given in Table 1. So over 95% of all the materials in the cloud is contained within a circle with a radius equal to 2.5. 2.3.3. Release conditions The continuous discharge of material is represented by the release of discrete clouds of material, and each cloud has an initial area equal to the area of the extended source. The interval between cloud releases is determined by the mean flow, and is equal to the time taken by the mean flow to transport the cloud a distance equal to the cloud diameter. So for a cloud with initial radius r0 in a mean flow U, the interval between cloud releases (TR ) is given by TR = r0 /|U|.
(16)
The initial quantity of material in each cloud (Q0 ) is then given by ˙ TR , Q0 = Q
(17)
˙ is the discharge rate. where Q 2.3.4. Motion of the cloud centroid After each cloud has been released into the turbulent velocity field its trajectory is calculated using a simple Euler integration, x(t + t) = x(t) + u(x, t)t, where the velocity at the cloud centroid (u(x, t)) is given by the sum of the KS velocity field and a random component u(x, t) = uks (x, t) + u (t − TR ) .
(18)
is a vector of random numbers drawn from a Gaussian distribution with zero mean and unit variance.
The random component models the effects of the small scale motions which lie outside the wavenumber range of the large scale simulation. The amplitude of the random velocity (u (t)) is the standard deviation of the cloud centroid velocities, calculated from an ensemble-average of the results from many smallscale simulations. The appropriate value for u depends on the elapsed time since the cloud was released, (t − TR ), rather than the absolute time. 2.3.5. Cloud growth and re-discretization All the clouds in the simulation are continually growing due to the action of small-scale motions. The rate at which the cloud radius increases depends on the time that has elapsed since the cloud was formed,
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
285
and is given by the results from KS1. Eventually, however, each cloud reaches the size of the grid cells in the concentration grid. At this point the fluid velocities in the large-scale simulation begin to vary on the same scale as the cloud itself, and contribute both to the spreading and the advection of the cloud. It is no longer reasonable to treat the cloud as a single entity, so all the clouds with radii equal to, or larger than, the grid spacing are split into smaller clouds. To avoid the creation of a large number of new clouds the existing clouds are actually re-discretized into a number of smaller clouds. The procedure for doing this is as follows: A grid is placed over the plume so as to cover all the clouds with radii greater than the splitting radius. The grid spacing is equal to the spacing of the concentration grid, although the splitting and concentration grids will not, in general, coincide. At the end of the splitting process one new cloud will be formed in each cell (except for those cells where the total quantity of material falls below some threshold—when the value is less than 10−10 Q, where Q is the quantity of material in the cloud, see Eq. (14), even smaller threshold value was also chosen for the simulation but the results were the same). In turn, the contribution of each cloud to each of the cells in the splitting grid is calculated. The concentrations are evaluated at fixed time intervals during the simulation, and output as data files at the end of the simulation. It is also possible to compute mean and fluctuating concentrations by taking the ensemble-average of many realizations.
3. Determination of parameters In this study, the CDM requires the following parameters: • • • • •
Energy spectrum for the turbulent velocity field E(k). Mean flow velocity U. Turbulence intensity u . Size of extended source. Grid spacing and grid size for concentration calculations.
The best approach would be to measure the energy spectrum in the flow that is to be modeled, but this can be expensive and difficult, if not impossible. Experimental measurements are made of the longitudinal one-dimensional (1D) spectrum function 11 (k1 ), which can be related to E(k). Since a spectrum may not be available, it is necessary to be able to synthesize some sort of model spectrum for use in these simulations, and to do this we have used a combination of theoretical arguments and available data. For the simulations of 2D homogeneous isotropic turbulence a certain form of energy spectrum is needed. Two-dimensional turbulence contains the two basic ingredients of randomness and convective nonlinearity (Hinze, 1975), and some of the statistical hypotheses which have been proposed for 3D turbulence should be applicable to 2D motion. Since the enstrophy (defined as half the mean-squared vorticity) tends to cascade toward large wavenumbers through non-linear interactions between different scales of motion in spatially homogeneous 2D turbulence, it has been proposed by Kraichnan (1967), Batchelor (1969) and Leith (1968) (referred to as KLB) that the enstrophy flux should be independent of k at high wavenumbers and equal to the rate of cascade of the mean-squared vorticity . Then a dimensional argument, similar to the one used in 3D turbulence for the Kolmogorov energy cascade (Kolmogorov, 1941a, b), assumes that E(k) is a function of and k only. This leads to the conclusion
286
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
that 2D turbulence simultaneously exhibits a direct cascade of enstrophy to large wavenumbers, up to a dissipation wavenumber k , and an inverse cascade of energy to small wavenumbers, down to wavenumber k = 0. In the limit of small viscosity, the inverse cascade is thought to proceed indefinitely in time to ever-larger scales, transferring virtually all of the energy input to wavenumber zero. The direct cascade is thought to come into balance with viscosity, transferring virtually all of the enstrophy input to k , where it will then be dissipated. In the long-time limit, a quasi-steady state is reached, in which two inertial ranges are established. According to the KLB theory, the energy range, which is only quasi-steady, scales as k −5/3 . The enstrophy range, which is in absolute equilibrium, should be as k −3 . I.e., from the enstrophy cascade concept, the energy spectrum is given by E(k) ∼ 2/3 k −3 .
(19)
For low wavenumbers, the spectrum is proportional to k 3 , due to resonant non-local interactions of two energetic wavenumbers (Basdevant et al., 1978). This leads to a spectrum of the form E(k) =
g1 k 3 (g22 + k 2 )3
,
(20)
where g1 and g2 are constants which depend on the flow. There are a number of ways of obtaining estimates for g1 and g2 . Perhaps the easiest is to fit (20) directly to appropriate measured spectra, but this will not always be possible, and then it is necessary to be able to estimate g1 and g2 by order-of-magnitude arguments. The best way of doing this will be through the use of estimates for the turbulence energy and the integral length scale. The total turbulence energy is simply the integral of the energy spectrum over wavenumber space (see A.10):
∞ 1 2 g1 E(k) dk = 2 . (21) u = 2 4g2 0 The integral length scale is more difficult to calculate; in 3D homogeneous isotropic turbulence the relationship between energy spectrum and the integral length scale is given by Batchelor (1953): ∞
∞ −1 k E(k) dk E(k) dk . (22) L= 0
0
In 2D homogeneous isotropic turbulence the relationship has a more complex form:
∞ ∞
∞ 2 R11 (r) dr u1 = J0 (kr)E(k) dk dr. L= 0
0
(23)
0
This result is derived in Appendix A. Substituting the suggested form for the energy spectrum (20) into (23) we obtain L = /4g2
(24)
from which g1 =
2 u21
4L2
(25)
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
287
and g2 =
4L
.
(26)
So (20) and (23) together are sufficient to determine g1 and g2 , given suitable estimates for the total turbulence energy and the integral length scale of the turbulence. In those cases where appropriate velocity measurements have been made, it will be possible to construct a typical spectrum. Usually, what is actually measured is u1 (t) at a fixed location in space, and this is then converted to u1 (x1 ) by means of Taylor’s “frozen eddy” hypothesis x1 = U1 t, where U1 is the component of the mean flow in the x1 direction. By taking the Fourier transform of this signal it is possible to obtain the 1D spectrum function (k1 ). Then the energy spectrum E(k) can be obtained from
∞ 3 d (27) E(k) = −k (k1 ) k12 − k 2 dk1 . 3 dk k 1 The derivation of this result is presented in Appendix Appendix B. In the simulations described in Section 4 we have assumed a spectrum of the form given in (20), and we have used it to model dispersion in two types of marine flow. The first of these (S1) is dispersion in a tidal current near a coast; there have been many field measurements in this type of flow (e.g. Grant et al., 1962; Bowden and Fairbairn, 1956; Bowden and Howe, 1964) and on the basis of those results we have assumed a mean flow of 1 m/s, a turbulence intensity of 10% and an integral length scale of about 250 m. This will give a ‘plume’ type dispersion, in which advection by the mean flow will be significant. The second simulation (S2) is an attempt to model the dispersion which might occur at a depth of 1000 m in the ocean. There have not been as many field measurements in this type of flow so our estimates for suitable parameter values will be less reliable. Ewart and Beninder (1981) measured the dispersion of a dye cloud in the ocean at a depth of 1000 m. They observed that at scales of 200 m or less, there is little evidence of severe distortion of the patch by eddy activity. We have reason to expect, therefore an increase in p somewhere between the scales of 200 m and 30 km, reflecting a change in the physical process predominantly responsible for horizontal mixing. (Their variable p was a measure of the rate at which the diffusion increased; d2 /dt ∝ p .) On the basis of this we have assumed an integral length scale of 1000 m. This immediately raises a potential problem in the field measurement of suitable spectra. For if the mean velocity is of the order of 10 cm/s and the integral length scale is of the order of 10 km, velocity measurements will have to be carried out for a period of at least 30 h in order to resolve these very large scales. This is not a trivial task! The alternative is to attempt to measure the spatial structure of the velocity fields directly, through measurements of the fluid velocities at many locations simultaneously. Data on deep ocean currents in the Pacific are also scarce. Near the coast of Japan the main current is the Kuroshio (Pickard and Emery, 1990) with mean speeds of between 75 and 250 cm/s. However, this is a surface current, which does not extend to a depth of 1000 m; there is some evidence of a deeper counter-current (Pickard and Emery, 1990) but it is not conclusive. Firing (1987) describes the mean zonal velocity structure of the Pacific. He identifies three main features below 400 m which appear to be permanent—eastward North and South Intermediate Counter currents with cores at about 700 m depth
288
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Table 2 Properties of the flow fields Simulation S1
S2
1.0 0.0 0.2 250 3.95 × 10−5 3.14 × 10−3 55
0.2 0.0 0.2 1000 9.87 × 10−8 7.85 × 10−4 55
U (m/s) V (m/s) u (m/s) L (m) g1 g2 LR (m)
Mean flow in x-direction Mean flow in y-direction Turbulence fluctuation Length scale Energy spectrum; g1 k 3 /(g22 + k 2 )3 Radius of extended source
and maximum mean speeds of about 15–20 cm/s, and the westward South Equatorial Intermediate Current with its core at about 900 m depth and maximum mean speed of about 15 cm/s. On the basis of these figures we have assumed a current velocity of 20 cm/s and a turbulence intensity of 100%. This will produce a simulation which will be dominated by the diffusive effects of the turbulence. The size of the extended source will depend on the small-scale model used for the jet mixing model; we have made an order-of-magnitude estimate for this, of 100 × 100 m2 . This gives a release radius of √ r = 100/ ≈ 55 m. The parameters used in these two simulations are given in Table 2. Of course all these parameters are capable of refinement; either by field measurements at the locations of interest or by a more thorough search of existing relevant field measurements. We have simply adopted reasonable values to demonstrate the application and the capabilities of the CDM.
4. Simulation results 4.1. Program parameters The size and resolution of the grid on which the concentrations are to be calculated depends on the resolution required. In both S1 and S2 we have computed concentrations on a 128 × 128 grid. In S1 the concentration grid covers an area of about 100 × 100 km2 (LS = 105 m) giving a grid spacing of about 780 × 780 m2 ; in S2 the concentration grid covers a much smaller area—25 × 25 km2 —giving a grid spacing of about 200 × 200 m2 . The input parameters for KS1 are shown in Table 3.The number of modes (Nk ) should be at least 50 to provide a reasonable discretization of the energy spectrum (Fung et al., 1992). The number of particles in each realization (Np ) can also be chosen freely, within the range of about 800–1600. The number of realizations (NR ) should be large enough to ensure that the ensemble average of the simulations converges; we have found that a minimum of about 50 realizations is required for this (Fung et al., 1992). It is important that only one cloud be released in each realization if statistics are to be calculated, because the ensemble-average must be taken over independent realizations. The total computing time is proportional to the number of realizations.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
289
Table 3 Input parameters for KS1 Simulation
LG (m) k1 (m−1 ) kN (m−1 ) Nk Np
r0 (m)
t (s)
T (s) NR
S1
S2
780.0 4.0 × 10−3 8.0 × 10−2 50 1600 1.0 50 55 1.15 × 104 100
200.0 4.0 × 10−3 8.0 × 10−2 50 1600 1.0 40 200 5.5 × 104 100
Grid size k1 = 2/LG kN = 2/LR No. of modes (minimum 50) No. of particles (minimum 900) Frequency factor Initial radius of pollutant Time step Duration No. of realization (minimum 30)
Table 4 Input parameters for KS2 Simulation
LG (m) rs (m) u (m/s) u (m/s) k1 (m−1 ) kN (m−1 ) Nk t (s) TR TS1 TS2 T (s) NR
S1
S2
780.0 390.0 1.0 0.2 6.0 × 10−5 4.0 × 10−3 50 110 t 32t 65t 7.7 × 104 50
200.0 100.0 0.2 0.2 2.5 × 10−4 4.0 × 10−2 50 150 4 t 15t 50t 6.0 × 104 392
Grid size Splitting radius 10–100% of mean flow k1 = 2/LS kN = 2/LG Minimum 50 TR = 2LR /|U| Obtained from KS1 Obtained from KS1 Duration No. of realizations
The total duration of the simulation (T) should be sufficient to ensure that the cloud radius exceeds LG . It will be necessary to continue beyond this point so that all the clouds to be included in the ensembleaverage exceed this radius. The input parameters for KS2 are shown in Table 4.The time interval between releases TR is obtained from (16); the time step t depends on two parameters—the release interval and the time step needed for an accurate integration of the cloud trajectories. As a rough rule, we estimate that the clouds should not 1 travel, on average, more than about 10 th of a grid square in a single timestep, from which we obtain (U + u )t LG /10
290
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
from this, t = min(TR , LG /10(U + u )).
(28)
The time until the cloud first splits (TS2 ) and the time interval between splitting (TS2 − TS1 ) are both obtained from KS1; this is described in Section 4.3. The number of realizations should be sufficient for the ensemble-average statistics to converge; in S1 we found that about 50 realizations were sufficient to obtain stationary statistics. In S2 we needed many more realizations; this is because the diffusive effects (which vary from realization to realization) are much stronger than the advection effects. Also, the concentrations were computed at a much finer scale in S2 than in S1, so the spatial variability due to diffusion becomes more important—if we were to smooth the concentrations in individual realizations we would obtain a faster convergence to the mean, but we would also lose some of the details of the fluctuations. 4.2. Simulation results—the velocity field The energy spectra used in these calculations are shown in Figs. 2a and b, together with the wavenumber ranges for KS1 and KS2. Fig. 3a shows a typical velocity field generated by KS with the full wavenumber range (i.e. k1 < k < kN ) whilst Figs. 3b and c show the corresponding fields for KS1 and KS2. To provide data for comparison with a standard Gaussian plume model we have calculated the diffusivity of the full velocity field by releasing a large number of particles and calculating their dispersion with time. Then the diffusion coefficients in the x- and y-directions are given by Kx =
Np 1 d 1 (xi − xR )2 2 dt Np
and
Ky =
i=1
Np 1 d 1 (yi − yR )2 , 2 dt Np
(29)
i=1
where Np is the number of particles and (xR , yR ) is the release point. The values are obtained from an ensemble average of 600 realizations. For S1 we obtain K ∼ 50 m2 /s, and for S2, K ∼ 75 m2 /s. 4.3. Simulation results—KS1 The two statistical processes which are needed from KS1 are the growth of the cloud radius (12) as a function of time, and the variance of the centroid velocity (11) as a function of time. Figs. 4a and b show the cloud radius for S1 and S2, respectively, together with curves which have been fitted to the data by the method of least squares. The general form of the curves is a quadratic 2 = a0 + a1 t + a2 t 2 ,
(30)
and the calculated values for a0 , a1 and a2 are given in Table 5.Figs. 5a and b show the variance of the centroid velocity, together with curves which have been fitted to the data by the method of least squares. The general form of the curve is a straight line 2u = b0 + b1 t
and the values of b0 and b1 are given in Table 5.
(31)
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
291
Fig. 3. A typical velocity field generated by kinematic simulations with different wavenumber ranges: (a) large scales k1 < k < kG , (b) small scales kG < k < kN and (c) full scales k1 < k < kN .
To check the effect of truncating the wavenumber range of KS1 we have also calculated 2 in a velocity field (S1) with a full wavenumber range (i.e. k1 < k < kN ). The results for the full and truncated wavenumber ranges are shown in Fig. 6, from which it is clear that in this simulation the longer length scales only begin to affect the ensemble-average cloud growth once the radius exceeds about 640 m. Since the wavenumber cut-off occurs at k = 4.0 × 10−3 the statistical results are reasonably insensitive to this truncation.
292
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Fig. 4. From KS1, the growth of the cloud radius is plotted as a function of time (symbols) together with the fitted curve (solid line): (a) for S1 and (b) for S2. Table 5 Coefficients for the fitted curves Simulation S1 a0 a1 a2 b0 b1
1.51 × 103 4.42 × 10−1 2.87 × 10−3 1.97 × 10−1 −6.78 × 10−6
S2 1.51 × 103 6.79 × 10−1 2.90 × 10−4 5.30 × 10−2 −7.79 × 10−7
4.4. Simulation results—KS2 In this section we discuss calculated concentrations and concentration fluctuations for both S1 and S2. To visualize dispersion in a single realization of the flow field we present concentration distributions at selected times for S1 (Figs. 7a and b) and S2 (Figs. 10a and b). Ensemble-averaged concentrations and concentration fluctuations were computed for S1 (50 realizations, Figs. 8a and b and 9a and b ) and for S2 (392 realizations, Figs. 11a and b and 12a and b). The apparent rectilinear outline of the clouds is due to the use of a lower concentration threshold in both the calculations and the plotting. It is sufficient to simply set this threshold below the lowest concentration of interest. 4.4.1. Simulation 1—S1 Near the source the dispersion of material is dominated by advection, because the mean flow is large compared with the turbulence, and because the area over which the concentration is being averaged is
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
293
Fig. 5. From KS1, the variance of the cloud centroid velocity is plotted as a function of time (symbols) together with the fitted curve (solid line): (a) for S1 and (b) for S2. The different symbols represent two different numerical techniques for tracking the growth and the motion of a cloud—a particle-tracking method (open circles) and a boundary-tracking method (closed circles). Basically, these two different numerical methods give the same statistics (Malik et al., 1991).
Fig. 6. Comparison of the growth of the cloud radius with time, with a full wavenumber range (symbols) and with a truncated wavenumber range (solid line) in KS1.
294
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Fig. 7. Plots of the concentration contours in one realization of S1 at different times: (a) t = 300t and (b) 600t.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Fig. 8. Plots of the mean concentration contours in S1 at different times: (a) t = 300t and (b) 600t.
295
296
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Fig. 9. Plots of the concentration fluctuation contours in S1 in different times: (a) t = 300t and (b) 600t.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Fig. 10. Plots of the concentration contours in one realization of S2 at different times: (a) t = 300t and (b) 600t.
297
298
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Fig. 11. Plots of the mean concentration contours in S2 at different times: (a) t = 300t and (b) 600t.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Fig. 12. Plots of the concentration fluctuation contours in S2 in different times: (a) t = 300t and (b) 600t.
299
300
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
large compared with the size of the plume. Consequently there is little evidence of the ‘meandering’ which would be observed in reality. As the time increases the concentration distribution in a single realization shows some “patchiness”—particularly within the plume—but the overall spread appears to be fairly smooth and symmetrical. However, if the instantaneous concentrations are compared with the corresponding ensemble-averaged concentrations (Figs. 8a and b) it is clear that there is quite a lot of internal structure to the plume. The plots of concentration fluctuations have been normalized on the local ensemble-averaged concentrations, and this emphasizes the fluctuations near the edges of the plume where the average concentrations are low. Within the plume the fluctuations are roughly symmetric about the centerline, although some large-scale variation is still visible. Even averaged over 50 realizations the outline of the plume shows some structure. In the S1 case, a higher resolution grid 256 × 256 was used to test the influences of the grid spacing on the computed mean concentration, the result was less than 5% difference in the maximum mean concentration value. Therefore, we can conclude that the current grid size is fine enough. 4.4.2. Simulation 2—S2 In S2 the turbulence intensity is equal to the mean velocity, so the material spreads out as a slightly elongated ‘cloud’ rather than as a plume. The spatial resolution of the concentration grid is much higher than that in S1, so it is possible to see much more of the structure of the cloud. Close to the source the plume meanders strongly, because most of the turbulent energy is at length scales which are large compared with the size of the cloud (Fig. 10a). At longer times the material begins to disperse more, creating several larger clouds (Fig. 10b); this gives rise to a very patchy concentration distribution. If the concentrations are averaged over a large number of realizations (Figs. 11a and b) a much smoother, larger cloud emerges. This demonstrates particularly well how an ensemble-averaged solution would obscure important details of the small-time, small-distance behavior of a single release, by predicting far too much dispersion close to the source. Even averaged over 392 realizations, however, the cloud still shows some structure. As with the results from S1 the highest fluctuations (when normalized on the local mean concentration) occur at the edges of the cloud. In this case, the plots of concentration fluctuations (Figs. 12a and b) emphasize the variability of the cloud from realization to realization. 4.5. Simulation results—a comparison between KS2 and a Gaussian plume As a simple check on the concentrations calculated by KS2 we have compared the ensemble-average concentration profiles in the x- and y-directions (for both S1 and S2) with those which would be obtained from a Gaussian plume in a uniform flow. The steady-state solution to the advection–diffusion equation (2) for continuous discharge into a uniform mean flow is given by ˙ Q (y − yR )2 C(x, y) = √ , (32) exp − 4(K/U )(x − xR ) 4(K/U )(x − xR ) ˙ into a flow with mean velocity U and diffusivity K. where material is released from (xR , yR ) at a rate Q Now in S1 a cloud with 1 ‘unit’ of material was released every 110 s, at (0, 0), into a flow with mean velocity 1 m/s and the computations of KS with the full wavenumber range give K ∼ 50 m2 /s. In S2,
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
301
1 unit was released every 1600 s, into a flow with a mean velocity of 0.2 m/s and diffusivity of 75 m2 /s. Substituting these into (32) we can compute the concentration along the centerline (i.e. y = 0) and the concentration profiles at different locations (x = constant). The steady-state Gaussian plume solution and the ensemble-average KS2 results are plotted together in Figs. 13 (S1) and 14 (S2). The agreement between the two is very good, except close to the source (where the Gaussian plume is unrealistic) and at the very end of the plume, where the dispersion has not yet reached a steady state. This result demonstrates that the CDM can reproduce the mean dispersion of a cloud of material whilst generating individual realizations which differ significantly from the mean. 4.6. Comparison with observations The wavenumber spectrum of horizontal velocity in the ocean is not normally determinable, requiring for direct calculation, elaborate current meter arrays. However, Zang and Wunsch (2001) based on the current meters in the North Pacific and estimated the zonal and meridional component of horizontal velocity, based on the estimated function (k, , , n, , ), where (k, ) are local horizontal wavenumbers; is frequency; (, ) are latitude and longitude and n is the vortical mode of variability. They found that the zonal-wavenumber spectrum of zonal velocity is proportional to k −1/2 at wavelength longer than 400 km, k −5/2 at wavelength between 100 and 400 km and k −4 at wavelength shorter than 100 km. We use this field measured spectrum for our KS simulation in order to compare with the tracer-release experiment data sets published by Sundermeyer and Ledwell (2001). In their work, tracer-release studies from coastal mixing and optics (CMO) dye experiments were used to examine the rates and mechanisms of lateral dispersion in the coastal ocean (which is in a stratified environment). Dispersion on scales of a few hundred meters to 10 km and for periods ranging from 2.5 to 5 days was investigated. Reversible and irreversible components of the lateral dispersion of the dye were defined by Sundermeyer and Ledwell (2001), using a simple model of vertical shear and horizontal strain. They found that the total diffusivities (reversible plus irreversible) estimated from the vertical integral of tracer ranged from K = 2.5 to 12.7 m2 /s. Irreversible diffusivities estimated from the dye along the target density surface were somewhat smaller and ranged from Kirrev = 0.3 to 4.9 m2 /s. In light of these results, we would like to see whether our KS simulation result could estimate a similar order of magnitude of the horizontal diffusivity. Parameters for this KS are k1 = 3 × 10−4 , kN = 0.1 and kG = 0.01 and LG = 2 km and similar to the works in, we found that the horizontal diffusivity is approximately equal to 5.4 m2 /s which is consistent with the results of Sundermeyer and Ledwell (2001).
5. Conclusion We have presented here a simulation of an unsteady random velocity field, whose statistics are consistent with the experimentally measured values of two-point space-time Eulerian and Lagrangian statistics of homogeneous turbulence (Fung et al., 1992). This KS should be useful in many studies which require such a field, and where there is no model for the small scale at present—it is a kind of “small-eddy simulation” that might complement large-eddy simulation. The simulation is not equivalent to that of a turbulent flow in two important respects: the vortical regions are not sufficiently elongated and the higher-order statistics are too closely Gaussian (Fung et al., 1992). The reason is that the simulation does not represent
302
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Fig. 13. (a) The steady-state Gaussian plume solution (thin line) along the centerline (Ny = 64) plotted together with the concentration calculated by KS2 (thick line) after 50 realizations of S1. (b) The Gaussian plume concentration profiles at different locations of x (thin line) plotted together with the concentration calculated by KS2 (thick line) after 50 realizations of S1. (i) x = 10.8 km, y = 1.15 km and centerline concentration Cm = 3.26, (ii) x = 26.5 km, y = 1.90 km and Cm = 1.85, (iii) x = 42.2 km, y = 2.41 km and Cm = 1.66, (iv) x = 58.0 km, y = 3.22 km and Cm = 1.12, (v) x = 73.7 km, y = 3.86 km and Cm = 8.34 × 10−1 , (vi) x = 89.4 km, y = 6.60 km and Cm = 2.73 × 10−5 , respectively. (c) The spread y of a Gaussian plume at different locations of x is plotted (thin line) together with the spread y calculated by KS2 after 50 realizations of S1.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
303
Fig. 14. (a) The steady-state Gaussian plume solution (thin line) along the centerline (Ny = 64) plotted together with the concentration calculated by KS2 (thick line) after 50 realizations of S2. (b) The Gaussian plume concentration profiles at different locations of x (thin line) plotted together with the concentration calculated by KS2 (thick line) after 50 realizations of S2. (i) x = 2.73 km, y = 2.67 km and Cm = 0.065, (ii) x = 6.76 km, y = 3.12 km and Cm = 0.054, (iii) x = 10.8 km, y = 3.90 km and Cm = 0.029, (iv) x = 14.8 km, y = 4.96 km and Cm = 0.007, respectively. (c) The spread y of a Gaussian plume at different locations of x is plotted (thin line) together with the spread y calculated by KS2 after 50 realizations of S2.
304
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
the dynamical processes in turbulence, which for example stretch vortices and affect higher-order statistics by making the turbulence intermittent. But the fact that the Lagrangian and pressure statistics agree with measurements suggested that these detailed aspects of the flow may not be so important when modeling many aspects of dispersion, combustion, two-phase flow, etc. Mean and fluctuating concentrations of pollutants has been computed using a new technique called the CDM. This model has been applied to simulate the dispersion of an uniform flow and the results are compared well with analytic, steady state solution. Furthermore, a preliminary attempt has been made to compare the rate of dispersion predicted by the model with the field data from the CMO dye studies and the computed diffusivity by the CDM has the same order of magnitude as that of the field experiments. This conclusion lends credence to the CDM. There are several ways in which the current model could be improved and developed. The most obvious weakness is the energy spectrum. We have derived a form for this, based on theoretical arguments, and we have estimated the values of the unknown coefficients. We have also applied a measured energy spectrum to our model but we anticipate that ultimately regional spectral models to describe different parts of the ocean around the world will emerge (Zang and Wunsch, 2001). At the moment the clouds in the CDM are all assumed to be circular, on the basis that they are modeling the average effect of the small scales, and that these produce, on average, radially symmetric concentration distributions. But in any one simulation the clouds will be far from circular, and this will lead to steepening of the concentration gradients. If it was felt that this effect might be important, it would be possible to develop a version of the CDM with non-circular clouds, using information from KS1 for higher moments of the cloud shape, and information from KS2 on local shear to suggest the direction in which elongation should occur. At the moment the CDM only simulates the dilution of the discharge, but because it computes individual realizations of a flow, it would be possible to include chemical reactions and two-phase behavior and these are the subjects of ongoing investigation. Acknowledgements Support from the Hong Kong Research Grant Council, China (Project nos. N/HKUST630/04 and 601203) and HIA02/03.SC04 are gratefully acknowledged. Appendix A. Let u(x) be a homogeneous random velocity field with a correlation function R(r). R(r) can be expressed by a Fourier transform over space, i.e.
Rij (r) = ij ()e−.r d, (A.1) where d is a 3D volume element at the position in wavenumber space and the integration is over all wavenumber space. The function ij () is called the spectral density. Since for real fields, Rij (r) = Rij (−r), then ij () = ij (−). In this case, (A.1) can be simplified to
Rij (r) = ij () cos(.r) d. (A.2)
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
305
However, for a 2D isotropic velocity field we can rewrite (A.2) in polar coordinates:
∞
Rij (r) =
2
0
cos(kr cos ) d ij (k)k dk,
(A.3)
0
where the direction of the vector is selected as the polar axis, i.e. for isotropic field ij () = ij (k), where k = ||. Also Rii (r) = 2Rii (r) (Batchelor, 1953). Therefore, (A.3) can be written as
∞ Rii (r) = J0 (kr) ii (k)k dk, (A.4) 0
where J0 (x) is the Bessel function, whilst k = (k12 + k22 )1/2 and r = (r12 + r22 )1/2 . Inversely
1 ∞
ii (k) = J0 (kr)Rii (r)r dr,
(A.5)
0
which for isotropic turbulence is related to the wavenumber–magnitude energy spectrum function E(k) by
ij () =
E(k) 2 (k ij − ki kj ). k 3
(A.6)
Substituting (A.6) into (A.4), we obtain
∞ R11 (r) = J0 (kr)E(k) dk.
(A.7)
0
For our model spectrum of the form E(k) =
g1 k 3 (g22 + k 2 )3
,
(A.8)
we have
∞
R11 (r) = g1 0
= g1
J0 (kr)k 3 (g22
+ k 2 )3
∞
dk = g1 0
J0 (kr)k (g22
r r2 K1 (g2 r) − K2 (g2 r) , 2g2 8
+ k 2 )2
dk
− g1 g22
∞ 0
J0 (kr)k (g22 + k 2 )3
dk (A.9)
where K1 (x) and K2 (x) are modified Bessel functions. The above integrals are carried out by using the formula
∞ Jp (bx)x p+1 a p−q bq Kp−q (ab), dx = 2q (q + 1) (x 2 + a 2 )q+1 0 where −1 < Re(p) < Re(2q + 3/2), a > 0 and b > 0. Putting b = r, a = g2 , p = 0 and q = 1 for the first integral and b = r, a = g2 , p = 0 and q = 2 for the second integral we obtain (A.9).
306
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
If the contribution to the energy u2 from the part of the wavenumber space between circles of radii k and k + dk is expressed as E(k)dk, then the total kinetic energy per unit mass of fluid is
∞
∞ k3 g1 1 2 u = E(k) dk = g1 dk = 2 . (A.10) 3 2 2 2 4g2 (g2 + k ) 0 0 The above integral is carried out by using the formula /
∞ ( /)(1 + n − /) p 1 x −1 , dx = n+1 n+1 q (n + 1) p (p + qx ) 0 where 0 < / < n + 1. Putting = 4, n = 2, = 2, p = g22 and q = 1 we obtain (A.10). From (A.9) and (A.10), we have the normalized longitudinal velocity correlation function f (r): f (r) = R11 (r)/u21 = 2g2 rK 1 (g2 r) + g22 r 2 K2 (g2 r)/2. Therefore, the turbulence integral length scale is given by
∞
∞ g22 ∞ 2 L= f (r) dr = 2g2 rK 1 (g2 r) dr + r K2 (g2 r) dr = . 2 0 4g2 0 0
(A.11)
(A.12)
The above integral is carried out using the formula
∞ 1++ 1+− −1 −−1 , x K (ax) dx = 2 a 2 2 0 where Re( + 1 + ± ) > 0, Re(a) > 0. Putting = 1, = 1 and a = g2 for the first integral and = 2, = 2 and a = g2 for the second integral we obtain (A.12). Appendix B. If we take a homogeneous and isotropic field u(x) on a straight line, then we will obtain the stationary random function of one variable u(x), for which the longitudinal 1D spectrum function (k1 ) may be defined as the Fourier transform of the auto-correlation function f (r), i.e.,
∞ 1 (k1 ) = u2 f (r1 ) cos(k1 r1 ) dr1 2 −∞
∞ =
11 () dk2 −∞
∞ E(k) =2 (k 2 − k12 ) dk2 . (B.1) k 3 0 Since k22 = k 2 − k12 , we have k2 dk2 = k dk (for fixed k1 ), therefore
2 ∞ E(k) 2 (k1 ) = (k − k12 )1/2 dk. k1 k2
(B.2)
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
307
Differentiation of (B.2) with respect to k1 gives
d(k1 ) 2 ∞ E(k) 1 2 (k − k12 )−1/2 (−2k1 ) dk, = dk1 k1 k2 2 therefore, we have 2 1 d(k1 ) =− k1 dk1
∞ E(k) k1
k
1 dk. k (k 2 − k12 )
(B.3)
Note that the solution of the integral equation (Uberoi and Kovásznay, 1953)
∞ () (x) = d x 2 − x 2 is
∞
d x [x (x)] · dx dx x 2 − 2
2 ∞ d2 [x (x)] · x 2 − 2 dx. = dx 2
() = −
2
(B.4)
Using (B.4), we can obtain the solution of (B.3) as follows:
2 E(k) 2 ∞ d2 1 d(k1 ) − = k · k12 − k 2 dk1 . 1 k k dk12 k1 dk1 Hence the corresponding explicit expression for E(k) in terms of the measurable function (k1 ) is
∞ 3 d (k1 ) · k12 − k 2 dk1 . (B.5) E(k) = −k dk13 k References Basdevant, C., Lesieur, M., Sadourny, R., 1978. Subgrid-scale modelling of enstrophy transfer in two dimensional turbulence. J. Atmos. Sci. 35, 1028–1042. Batchelor, G.K., 1949. Diffusion in a field of homogeneous turbulence I. Eulerian analysis. Aust. J. Sci. Res. 2, 437–442. Batchelor, G.K., 1952. Diffusion in a field of homogeneous turbulence II. The relative motion of particles. Proc. Cambridge Philos. Soc. 48, 345–362. Batchelor, G.K., 1953. The Theory of Homogeneous Turbulence. Cambridge University Press, Cambridge. Batchelor, G.K., 1969. Computation of the energy spectrum in homogeneous two-dimensional turbulence. Phys. Fluids Suppl. II 12, 233–239. Bowden, K.F., Fairbairn, L.A., 1956. Measurements of turbulent fluctuations and Reynolds stresses in a tidal current. Proc. R. Soc. London Ser. A 237, 422–438. Bowden, K.F., Howe, M.R., 1964. Observations of turbulence in a tidal current. J. Fluid Mech. 17, 271–284. Drummond, I.T., Münch, W., 1991. Distortion of line and surface elements in model turbulent flows. J. Fluid Mech. 225, 529–544. Drummond, I.T., Duane, S., Horgan, R.R., 1984. Scalar diffusion in simulated helical turbulence with molecular diffusivity. J. Fluid Mech. 138, 75–91.
308
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
Durbin, P.A., 1980. A stochastic model of two-particle dispersion and concentration fluctuations in homogeneous turbulence. J. Fluid Mech. 100, 279–302. Ewart, T.E., Beninder, W.P., 1981. An observation of the horizontal and vertical diffusion of a passive tracer in the deep ocean. J. Geophys. Res. 86, 10974–10982. Firing, E., 1987. Deep zonal currents in the central equatorial Pacific. J. Marine Res. 45, 791–812. Flohr, P., Vassilicos, J.C., 2000. A scalar subgrid model with flow structure for large-eddy simulations of scalar variances. J. Fluid Mech. 407, 315–349. Fung, J.C.H., 1993. Gravitational settling of particles and bubbles in homogeneous turbulence. J. Geophys. Res.—Oceans 98, 20,287–20,297. Fung, J.C.H., 1998. Effect of nonlinear drag on the settling velocity of particles in homogeneous isotropic turbulence. J. Geophys. Res.—Oceans 103, 27,905–27,917. Fung, J.C.H., Vassilicos, J.C., 1998. Two-particle dispersion in turbulent-like flows. Phys. Rev. E 57, 1677–1690. Fung, J.C.H., Hunt, J.C.R., Perkins, R.J., Wray, A.A., Stretch, D., 1991. Defining the zonal structure of turbulence using the pressure and invariants of the deformation tensor. In: Johansson, A.V., Alfredsson, P.H. (Eds.), Proceedings of the Third European Turbulence Conference, Stockholm, 1990. Springer, Berlin. Fung, J.C.H., Hunt, J.C.R., Malik, N.A., Perkins, R.J., 1992. Kinematic simulation of homogeneous turbulent flows generated by unsteady random Fourier modes. J. Fluid Mech. 236, 281–318. Fung, J.C.H., Hunt, J.C.R., Perkins, R.J., 2003. Diffusivities and velocity spectra of small inertial particles in turbulent-like flows. Proc. R. Soc. London Ser. A 459, 445–493. Garrett, C., 1983. On the initial streakiness of a dispersing tracer in two- and three-dimensional turbulence. Dyn. Atmos. Oceans 7, 265–277. Gotoh, T., 2002. Small-scale statistics of turbulence at high Reynolds numbers by massive computation. Comput. Phys. Commun. 147, 530–532. Grant, H.L., Stewart, R.W., Moillet, A., 1962. Turbulent spectra from a tidal channel. J. Fluid Mech. 12, 241–268. Hinze, J.O., 1975. Turbulence. McGraw-Hill, New York. Kolmogorov, A.N., 1941a. The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Proc. Acad. Sci. USSR 30, 301–305. Kolmogorov, A.N., 1941b. On degeneration of isotropic turbulence in an incompressible viscous liquid. Proc. Acad. Sci. USSR 31, 538. Kraichnan, R.H., 1967. Inertial range in two-dimensional turbulence. Phys. Fluids 8, 575–598. Kraichnan, R.H., 1970. Diffusion by a random velocity. Phys. Fluids 13, 22–31. Leith, C.E., 1968. Diffusion approximation for two-dimensional turbulence. Phys. Fluids 11, 671–673. Luhar, A.K., Britter, R.E., 1989. A random walk model for dispersion in inhomogeneous turbulence in a convective boundary layer. Atmos. Environ. 23 (9), 1911–1924. Malik, N.A., 1990. Distortion of material surfaces in steady and unsteady flows. In: Proceedings of the IUTAM Symposium on Topological Fluid Mechanics, CUP, pp. 617–627. Malik, N.A., Vassilicos, J.C., 1999. A Lagrangian model for turbulent dispersion with turbulent-like flow structure: comparison with direct numerical simulation for two-particle statistics. Phys. Fluids 11, 1572–1580. Malik, N.A., Perkins, R.J., Hunt, J.C.R., 1991. Computational modelling of dispersion in coastal, estuarine and river regions with application to sediment transport. Report to WRC. Maxey, M.R., 1987. The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields. J. Fluid Mech. 174, 441–465. Nicolleau, F., Vassilicos, J.C., 2003. Turbulent pair diffusion. Phys. Rev. Lett. 90, 024503. Pickard, G.L., Emery, W.T., 1990. Descriptive Physical Oceanography. fifth ed. Pergamon Press, Oxford. Rogallo, N.M., Moin, P., 1984. Numerical simulation of turbulent flows. Annu. Rev. Fluid Mech. 16, 99–157. Sawford, B.L., 1985. Lagrangian statistical simulation of concentration mean and fluctuation fields. J. Clim. Appl. Meteorol. 24, 1152–1166. Shraiman, B.I., Siggia, E.D., 2000. Scalar turbulence. Nature 405, 639–646. Sundermeyer, M.A., Ledwell, J.R., 2001. Lateral dispersion over the continental shelf: analysis of dye release experiments. J. Geophys. Res. 106, 9603–9621. Thomson, D.J., 1984. Random walk modelling of diffusion in inhomogeneous turbulence. Q. J. R. Meteorol. Soc. 110, 1107–1120.
J.C.H. Fung, R.J. Perkins / Fluid Dynamics Research 40 (2008) 273 – 309
309
Thomson, D.J., 1986. A random walk model of dispersion in turbulent flows and its application to dispersion in a valley. Q. J. R. Meteorol. Soc. 112, 511–530. Thomson, D.J., 1987. Criteria for the selection of the stochastic models of particle trajectories in turbulent flows. J. Fluid Mech. 180, 529–556. Thomson, D.J., 1990. A stochastic model for the motion of particle pairs in isotropic high Reynolds number turbulence, and its application to the problem of concentration variance. J. Fluid Mech. 210, 113–153. Turfus, C., Hunt, J.C.R., 1987. A stochastic analysis of the displacements of fluid elements in inhomogeneous turbulence using Kraichnan’s method of random modes. In: Comte-Bellet, G., Mathieu, J. (Eds.), Advances in Turbulence. Springer, Berlin, pp. 191–203. Uberoi, M.S., Kovásznay, L.S.G., 1953. On mapping and measurement of random fields. Quart. Appl. Math. 10, 375–393. van Dop, H., Nieuwstadt, F.T.M., Hunt, J.C.R., 1985. Random walk models for particle displacements in inhomogeneous unsteady turbulent flows. Phys. Fluids 28 (6), 1639–1653. Vassilicos, J.C., Fung, J.C.H., 1995. The self-similar topology of passive interfaces advected by two-dimensional turbulent-like flows. Phys. Fluids 7, 1970–1998. Vincent, V., Meneguzzi, M., 1991. The spatial structure and statistical properties of homogeneous turbulence. J. Fluid Mech. 225, 1–20. Yeung, P.K., 2002. Lagrangian investigations of turbulence. Annu. Rev. Fluid Mech. 34, 115–142. Zang, X., Wunsch, C., 2001. Spectral description of low-frequency oceanic variability. J. Phys. Oceanogr. 31, 3073–3095.