A&A 421, 775-782 (2004)
DOI: 10.1051/0004-6361:20035606
B. Dintrans^{1} - A. Brandenburg^{2}
1 - Observatoire Midi-Pyrénées, CNRS UMR5572, 14 avenue
Edouard Belin, 31400 Toulouse, France
2 - NORDITA,
Blegdamsvej 17, 2100 Copenhagen Ø, Denmark
Received 31 October 2003 / Accepted 2 April 2004
Abstract
The excitation of internal gravity waves by an entropy bubble
oscillating in an isothermal atmosphere is investigated using direct
two-dimensional
numerical simulations. The oscillation field is measured by a
projection of the simulated velocity field onto
the anelastic solutions of the linear eigenvalue problem for the
perturbations. This
facilitates a quantitative study of both the spectrum and the amplitudes
of excited g-modes.
Key words: hydrodynamics - waves - stars: oscillations - stellar dynamics
The problem of the excitation of internal gravity waves (hereafter IGWs) in stably stratified media is often studied in connection with the possible detection of solar g-modes (see e.g. the latest attempt of g-mode detection in the GOLF data by Gabriel et al. 2002). As an example, two- and three-dimensional numerical simulations of penetrative convection have shown that it is possible to excite such waves below the convection zone from the penetration of strong downward plumes into the stable radiative zone located below (Hurlburt et al. 1986; Hurlburt et al. 1994; Brandenburg et al. 1996; Brummell et al. 2002; Dintrans et al. 2003).
Furthermore, IGWs play a major role in stellar evolution as they can transport angular momentum and/or chemical elements over long distances in the stellar interior. This transport mechanism has been proposed to explain the quasi-solid rotation profile of the solar core as revealed by helioseismology (Kumar et al. 1999) as well as the lithium depletion observed in low mass stars (Talon & Charbonnel 1998). However, the correct amount of IGWs excited by, say, penetrative convection remains still unknown as numerical studies were essentially led from a qualitative point of view, so that the mode amplitudes were not deduced.
The main goal of this paper is to show that it is possible to determine quantitatively the amplitudes of gravity waves propagating in a stable zone of a numerical simulation using the anelastic subspace. This subspace is built from the solutions of the associated anelastic linear eigenvalue problem for the perturbations. The projection of the simulated velocity field onto this basis yields the IGW amplitudes. We will present here the application of this method to the simple problem of g-mode oscillations in an isothermal atmosphere. The dynamics of such an atmosphere is well known (Brandenburg 1988) so that this problem is ideal to illustrate and test the validity of our method, before applying it to the more difficult problem of IGWs excited by penetrative convection, where the atmosphere is in general non-isothermal.
After presenting the hydrodynamical model describing our isothermal atmosphere (Sect. 2), we apply two classical and widely used methods to detect the g-modes excited by an oscillating entropy bubble and show their limitations (Sect. 3). Next, we introduce our new method based on the anelastic subspace and give the analytic forms we found for both the eigenfrequencies and eigenvectors of the anelastic set of linear equations for the perturbations (Sect. 4). We then apply this formalism to the same simulation used to test the two classical methods discussed in Sect. 3 and show that we now have access to both the spectrum and amplitudes of the IGWs (Sect. 5). In particular, we illustrate the usefulness of our method by studying the influence of the box geometry on the mode amplitudes. Finally, we conclude in Sect. 6 by introducing the next step of this work, which will consist of the application of the anelastic subspace method to the detection of IGWs that are stochastically excited by penetrating convective plumes.
We consider the two-dimensional propagation of internal gravity waves
in an isothermal atmosphere consisting of a layer of depth dof a perfect gas at temperature T_{0} embedded between two horizontal
impenetrable boundaries. The fluid is non-rotating and stratified with
constant gravity and its properties like its kinematic viscosity, ,
and specific heats, c_{p} and c_{v}, are assumed to be constant (with
for a monatomic gas). Initially, the layer is in
hydrostatic equilibrium such that its density is given by
We choose the depth of the layer d as the length unit,
as
the time unit
and
and T_{0} as the units of density and temperature, respectively. The evolution
of this layer is governed by the equations of conservation of mass,
momentum and energy
When imposing the boundary conditions, we assume the top and bottom
surfaces to be stress-free with the fixed temperature (or,
equivalently, internal energy), that is
The fully nonlinear set of Eqs. (2) is solved using the hydrodynamical code described in Nordlund & Stein (1990) and Brandenburg et al. (1996). The time advance is explicit and uses a third order Hyman scheme (Hyman 1979). All spatial derivatives are calculated using sixth order compact finite differences (Lele 1992). The two-dimensional simulations presented here were done with the same resolution (i.e. 256 zones in the horizontal direction and 300 in the vertical) and aspect ratios, , between 2 and 6. Here, L_{x} and L_{z} denote the domain size in the horizontal and vertical directions, respectively. Moreover, the dimensionless sound speed and the gravitational acceleration are also the same for all runs and set equal to . As a consequence, the dimensionless pressure scale height is H=3/5 (i.e. the box extends over 1.66H in the vertical direction so that the density contrast between bottom and top is around ), while the Brunt-Väisälä frequency is . Finally, the dimensionless kinematic viscosity has been set equal to in all simulations, except in the study of the influence of the viscosity on the mode damping rates (Sect. 5.2) where we also considered smaller values down to .
As initial condition to excite IGWs in the layer, we choose the
perturbation to be an entropy bubble (Brandenburg 1988). For a given
entropy perturbation, it is indeed easy to find the approximate position
of the point around which the bubble will begin to oscillate
and thus generate IGWs. In dimensionless units we have
Figure 1: Velocity field superimposed on a grey scale representation of the entropy perturbation for a two-dimensional simulation of an entropy bubble oscillating in an isothermal atmosphere. The asterisks shown in panel a) are used in Fig. 2a. | |
Open with DEXTER |
Figure 1 shows an example of such a two-dimensional hydrodynamical simulation where the velocity field has been superimposed onto a grey scale representation of the entropy. At t=0, the negative entropy perturbation looks like a bubble at x=0and z=0.2 (Fig. 1a) and no velocity field is present. Under the effect of gravity, the bubble begins to descend (Fig. 1b) and stabilizes at the depth z=0.5. The bubble thus oscillates around this equilibrium position and IGWs are generated in the whole domain (Fig. 1c). The question is what is the spectrum and amplitude of the internal waves excited by this oscillating bubble.
Two techniques are commonly used to measure wave fields propagating in direct hydrodynamical simulations:
Figure 2: Detection of IGWs using method M1. a): time evolution of the vertical velocity recorded at three different positions in the domain (denoted by stars in Fig. 1a). b): corresponding temporal power spectra. Peaks below the buoyancy frequency N=0.82 correspond to internal gravity waves whereas the peak around is the fundamental acoustic mode. | |
Open with DEXTER |
Figure 3: Detection of IGWs using method M2. Left: temporal power spectra in the -plane for three different degrees . Right: the resulting spectra after an integration in depth. Dotted lines shown in all plots correspond to the bounding frequency , above which g-modes cannot exist. | |
Open with DEXTER |
In Fig. 3 we summarize the application of the second method M2 to detect IGWs in the same simulation of Fig. 2. We now have access to more informations than with method M1. At first, the diagram confirms: (i) that the mode around is indeed a radial one and that it corresponds to the fundamental acoustic mode as no radial node is visible; (ii) that gravity modes are strictly nonradial modes as no peaks are visible below the dotted line .
Second, gravity modes with similar frequency but different degrees are now well separated. As an example, the fundamental g-mode at and the first overtone at have almost the same frequency, . Using the -plane for different 's allows us to separate modes and to show that one mode corresponds to a fundamental one while the other one is a first overtone (note the two "bumps'' around in the -diagram for ).
The classical methods M1 and M2 presented above have some disadvantages which render them inappropriate for a careful detection of IGWs in hydrodynamical simulations:
Our new method, which eliminates these shortcomings, is inspired by
the work of Bogdan et al. (1993, hereafter BCM) who have
developed a tool to measure the acoustic emission generated in
their simulations of turbulent convection. They extracted the
acoustic field by projecting the convection field onto the acoustic
subspace build from the eigenfunctions of the associated linear
oscillations problem. Here we are interested in gravity waves rather
than sound waves and this allows us to adopt a simpler procedure than
theirs. Thus, we project our simulation data onto the anelastic
subspace, that is, we filter out the acoustic waves in the linear system
of oscillation equations. For a time dependence of normal modes of the
form
,
and in the ideal limit ,
the anelastic
set of equations reduces to (Dintrans & Rieutord 2001; Rieutord &
Dintrans 2002)
Because we adopt periodic boundary conditions in the x-direction,
we can seek solutions of the form
As shown in Appendix A, the eigenfunctions
which are
solutions of the eigenvalue problem (5) are orthogonal,
that is we have
Figure 4: Vertical structure of the fundamental anelastic mode (dark line) and first overtone (grey line) at (a) and (b). Solid lines mark the horizontal displacement while dot-dashed lines are for the vertical one . Eigenfunctions have been normalized by imposing . | |
Open with DEXTER |
Figure 4 shows four anelastic modes at (a) and (b) with n=0 (dark lines) and n=1 (grey lines). We remark that anelastic eigenfunctions hardly change with the -values as does not explicitly depend on k_{x}, while only the -amplitude involves k_{x}. As a consequence, it is the ratio which makes sense so that becomes more and more important compared to with increasing k_{x}, as observed in Fig. 4 when comparing solid to dot-dashed lines for the same order n. In the case of method M2, the use of the vertical mass flux (instead of, say, ) is justified a posteriori as a good proxy for detecting IGWs propagating in the simulation (see Fig. 3).
Figure 5: Upper: same as Fig. 3 at , except that amplitudes have been restored using third direction. Lower: comparisons between the normalized vertical mass flux profiles deduced from the upper power spectra (solid lines) with the ones obtained from the anelastic modes shown in Fig. 4a (dotted lines). | |
Open with DEXTER |
Once the anelastic eigenfunctions are known, one can compare the corresponding vertical mass flux with the ones deduced from the power spectra used in method M2. Indeed, the peaks appearing at a given frequency in Fig. 3 correspond to g-modes so that their vertical profiles can be related to the eigenfunctions of the linear modes. This was done by BCM for the case of acoustic modes propagating in their convection simulations.
The upper panel of Fig. 5 shows the depth dependence of the two peaks appearing at in Fig. 3, which correspond to the n=0 ( ) and n=1 ( ) modes, respectively. We thus recover what BCM called the "shark fin'' profiles, that is, the power at any depth is located around the theoretical anelastic frequency. Taking a mean profile around eigenfrequencies and allows us to compare the vertical mass flux observed in the simulation with those calculated from the linear anelastic modes. The lower panel of Fig. 5 shows that the agreement between the two profiles is quite remarkable, meaning that the anelastic modes reproduce well the long-period dynamics of the oscillating isothermal atmosphere.
Figure 6: Left: real parts of the complex mode amplitude for the three modes shown in Fig. 4. Right: corresponding temporal power spectra. Dotted lines mark the location of the anelastic eigenmodes and , respectively. Amplitudes on each plot have been multiplied by one thousand. | |
Open with DEXTER |
The left hand panel of Fig. 6 shows the real parts of the complex amplitudes and c_{21}, that is, the result of the projection of the simulation shown in Fig. 1 onto the three anelastic modes plotted in Fig. 4. By taking a Fourier transform of these sequences (Fig. 6, right), one obtains three peaks which agree remarkably well with the theoretical anelastic frequency given by (9), i.e. the complex amplitudes behave as
Figure 7: Time-evolution of the amplitude modulus for the three modes shown in Fig. 4. Grey lines superimposed correspond to an exponential decay of the form . | |
Open with DEXTER |
As we have access to the real and imaginary parts of , we now compute the amplitude modulus ; Fig. 7 shows its time evolution for the three modes in Fig. 6. A regression analysis of the curves allows us to determine the shape of the envelope and we found that follows a simple exponential decay law of the form . Because of the periodic behaviour (10), it means that the mode amplitude behaves like
Figure 8: Dependence of the time-averaged kinetic energy on degree and radial order n. | |
Open with DEXTER |
We will now show that one of the main advantages of the anelastic subspace resides in the fact that it allows to calculate the contribution of every g-mode to the time-averaged kinetic energy. We demonstrate in Appendix D that a Parseval-Bessel relation also exists in the anelastic subspace as
This relation illustrates the usefulness of the anelastic subspace when one wants to perform a detailed study of the IGW excitation. By using the classical Parseval-Bessel relation applied in Fourier space, it is already possible to find the amount of kinetic energy contained at a given wavenumber k_{x}, since
The main advantage of our anelastic subspace method is that it lifts this degeneracy by isolating contributions that come from g-modes only. Comparing to its classical form (13), the anelastic Parseval-Bessel relation (12) introduces the radial order n so that one can extract the contribution of the g-mode in the kinetic energy and then deduce which modes contribute most to it, meaning the modes that are the most excited by the oscillating bubble. This is what we did in Fig. 8 where we also plotted, for each , the contributions coming from the n=0, n=1 and n=2 g-modes. It allows us to show that the contribution is in fact almost entirely composed of the (1,0) and (1,1) g-modes, as they respectively contain around 33% and 15% of the time-averaged kinetic energy of this simulation. This kind of information is particularly relevant when one deals with the problem of wave excitation.
In order to assess the restriction imposed by the finite size of the box, we now study the influence of the box geometry onto the wave excitation. We therefore perform three runs with three different aspect ratios (2, 4 and 6) and calculate the g-mode contributions to the time-averaged kinetic energy in each case.
Figure 9 shows the spectral power in the plane for the three different aspect ratios. One recovers clearly the well-known dispersion relation for g-modes (e.g. Stein & Leibacher 1974); see also its anelastic counterpart in Eq. (9). One also verifies that increasing the horizontal extent of the box allows the power to be distributed among a larger number of horizontal modes. This phenomenon has also been observed in the problem of p-modes excitation by turbulent convection (Stein & Nordlund 2001; Nordlund & Stein 2001). Unfortunately, such diagrams are generally not easy to interpret when the system is not isothermal and the modes stochastically driven; see for example the three-dimensional simulations of Brandenburg et al. (1996) where the resulting turned out to be very noisy.
Figure 9: diagrams which three different aspect ratios (2, 4 and 6). Crosses correspond to the anelastic frequencies given by (9). | |
Open with DEXTER |
Contrasting this with the anelastic subspace method, it is straightforward to extract at a given wavenumber k_{x} the energy from modes with different radial order n. As expected, the spectral distribution of mode energies integrated over all n collapses onto an universal curve (Fig. 10). It turns out that the mode energy scales as for small enough values of k_{x}, i.e. for . However, both the exponent and the cut-off are not unique and depend on the size of the initial entropy bubble, as verified by running several simulations with bubble radii ranging from R=0.1 to R=0.4. We found that the larger the bubble, the smaller is the cutoff value of k_{x}, meaning that big bubbles preferentially generate g-modes at large scales. It would be interesting to investigate in more detail the origin of such behaviour, but this is clearly beyond the scope of this paper.
Figure 10: Distribution of the energy contained in g-modes for different aspect ratios (2, 4 and 6) in a log-log representation. | |
Open with DEXTER |
We have presented a new and accurate method to measure the internal gravity waves propagating in a numerical simulation. This is of prime importance when one deals with the problem of the generation of IGWs in a stably stratified medium such as the radiative zone of a star. The knowledge of the vigour of the velocity field in such a zone is the main input of any theory of wave transport of angular momentum and/or chemicals in stellar radiative zones (e.g. Schatzman 1993).
However, we have shown that classical methods based only on Fourier transforms in space and/or time of the simulated velocity field do not allow a quantitative determination of the mode amplitudes. The reason is mainly that the contribution of the propagating g-modes to the total energy in the simulation is diluted among the other ones due to p-modes and turbulent motions.
Our method is based on the projection of the simulated velocity field onto the anelastic subspace built from the solutions of the linear problem for the perturbations. Once the resulting time-dependent coefficients of the velocity field are computed, one can deduce the amplitudes of every g-mode and measure their kinetic energy at every timestep of the simulation.
We have applied this formalism to the measurement of IGWs excited by an entropy bubble oscillating in an isothermal atmosphere. We have shown that the mode amplitudes follow the same law than those of a linear damped harmonic oscillator, that is, the periodic amplitudes decrease exponentially with time, with a time scale that is simply proportional to the inverse of the viscosity; see Eq. (11). By using a Parseval-Bessel relation valid in the anelastic subspace, we have deduced the contribution of every g-mode to the time-averaged kinetic energy and have shown that large-scale motions are mainly composed of n=0 and n=1 modes. Such a result is in general impossible to obtain in a more complicated settings using a classical diagram.
The fact that the mode amplitudes in the anelastic basis follow the same law than those of a linear damped oscillator is important in the case of a stochastic excitation of IGWs, such as the one encountered when one deals with numerical simulations of penetrative convection. Indeed, in this case the mode amplitude evolves either chaotically around zero when the mode is not excited, or in a periodic fashion following Eq. (11) when the mode is excited. By computing time-frequency diagrams of that time-dependent amplitude, one can extract the time intervals during which the mode has been really excited. In a forthcoming paper, we will present the application of the anelastic subspace method we tested here to the problem of IGWs stochastically excited by penetrative convection (Dintrans et al. 2004).
Acknowledgements
This work has been supported by the European Commission under Marie-Curie grant No. HPMF-CT-1999-00411. Calculations were carried out on the CalMip machine of the "Centre Interuniversitaire de Toulouse'' (CICT) which is gratefully acknowledged. We also thank Michel Rieutord and the referee for useful comments.
We start from the complete set of equations written for two eigenmodes and (Dintrans & Rieutord 2001)
(A.1) |
(A.2) |
Finally, we integrate this last equation over the volume and use the Green-Ostrogradsky theorem to obtain
(A.3) |
(A.4) |
We start from the anelastic set of Eqs. (4)
(B.3) |
(B.5) |
(B.6) |
(B.8) |
Table B.1: Exact eigenfrequencies for the first five g-modes of degree , corresponding to an horizontal wavenumber (here L_{x}=2), in the complete case, in the anelastic approximation, and the corresponding relative errors.
The set of equations for the complete case is obtained by, first, linearizing the system (2) for infinitesimal adiabatic perturbations and, second, by assuming the time and horizontal dependences of the modes to be and or , respectively. We thus arrive at
Lamb (1909, 1932) found analytical solutions of the system (B.9) in the case of an isothermal atmosphere. By eliminating and in the system (B.9), we obtain a second-order differential equation for only,
Figure B.1: Comparison of the complete (solid lines) and anelastic (dashed lines) normalized eigenfunctions for and n=[0,1,2]. |
We summarize in Table B.1 the eigenfrequencies for the complete case deduced from (B.14) for the first five g-modes of degree and compare them with their anelastic counterparts given by the formulæ (B.7). The agreement is quite remarkable since the relative error done by the anelastic approximation is less than from the fundamental g-mode. Concerning the associated eigenvectors, we showed that the vertical displacements are the same in both cases so that only the horizontal ones differ. However, Fig. B.1 shows that the agreement on the normalized -eigenfunctions is also very good, except for the mode at n=0for which the condition (B.11) is not well fulfilled.
In conclusion, the filtering of acoustic waves from the dynamics of the isothermal atmosphere does almost not change the g-mode eigenfrequencies and eigenfunctions . The basis built from these anelastic eigenfunctions can therefore be used with good confidence to project out velocity fields from hydrodynamical simulations.
It should be pointed out that the application of the anelastic subspace method will only give good results when the acoustic and gravity spectra are well separated, that is, when the ratio is large. This is the case for the chosen isothermal setup, where the resulting anelastic eigenfrequencies/eigenvectors are very close to those of the complete problem (see Table B.1 and Fig. B.1). However, this ratio decreases with larger stratification, in which case one may be forced to solve numerically the oscillation Eqs. (B.9) of the complete problem instead of those of the anelastic subset (B.1), which does not however introduce any particular difficulty.
We show in this appendix that it is possible to find the relation
governing the damping rate of a mode and the viscosity using the work
integrals formalism. We start from the anelastic equations written for
the velocity
(C.1) |
We multiply the momentum equation by to obtain
(C.2) |
(C.3) |
(C.4) |
A similar relation also applies in the complete case, provided one includes the contribution of pressure fluctuations to the energy, that is,
In this appendix, we will show that it is possible to relate the total
kinetic energy embedded in the box to the one calculated in the anelastic
subspace. To do that, we need a similar relation than the Parseval-Bessel
one which states that the total energy is conserved when one works in
Fourier space, i.e.
(D.1) |
In our case, we only perform a 1-D Fourier transform of the velocity field in the horizontal direction (corresponding to a wavenumber ), so that the Parseval-Bessel relation for the total kinetic energy at a given time is
(D.3) |
= | |||
= | |||
(D.4) | |||
= | |||
= |
(D.5) |