Radial mixing in protoplanetary accretion disks

Radial mixing in protoplanetary accretion disks

VII. -dimensional transport of tracers
M. Wehrstedt Institut für Theoretische Astrophysik, Universität Heidelberg, Albert-Überle-Str. 2, 69120 Heidelberg, Germany (gail@ita.uni-heidelberg.de)    H.-P. Gail Institut für Theoretische Astrophysik, Universität Heidelberg, Albert-Überle-Str. 2, 69120 Heidelberg, Germany (gail@ita.uni-heidelberg.de)
Received XXX/ Accepted XXX
Key Words.:
Accretion, accretion disks – solar system: formation – dust, extinction
offprints: H.-P. Gail


Aims:The detection of significant concentrations of crystalline silicates in comets indicates an extensive radial mixing in the primordial solar nebula, i.e. the protoplanetary accretion disk of our solar system. In studying the radial transport of matter within protoplanetary disks by numerical model calculations it is essential to resolve the vertical disk structure since matter is mixed radially inward and outward by a complex -dimensional flow pattern within the disk that is superposed on the global inward directed accretion flow. It is further essential to follow numerically the advection-diffusion processes over a period of at least 10 yrs to allow for a full development of the radial concentration profile built up by radial mixing from the warm inner to the cool outer parts of protoplanetary disks beyond of 10 AU.

Methods:Numerical model calculations of protoplanetary accretion disks with radial and vertical mixing are performed by solving a set of -dimensional transport-diffusion-reaction equations for some important tracers self-consistently with the set of disk equations in the +-dimensional approximation. The global D velocity field of the disk is calculated from an approximate analytical solution for the meridional flow pattern, which exhibits an inward drift in the upper layers and an outward drift in the midplane in most parts of the disk. The disk model is based on the -prescription of viscosity and considers vertical self-gravitation of the disk. This kind of semianalytical approximations allows with presently available computer capacity to follow the evolution of the disk and the transport and mixing of tracers in vertically resolved 2D-models over the required long periods of disk evolution. The mixing processes in the disk are studied for the following species: amorphous silicate grains (forsterite, enstatite) which crystallise by annealing in the warm inner parts of the disk, and carbonaceous grains which are destroyed by surface reactions with molecules at elevated temperatures.

Results:Considerable fractions of crystallised silicates and methane (formed as a by-product of carbon combustion) are transported to the site of comet formation far from the protosun within a period of 10 yrs. The -dimensional transport of tracers in the solar nebula therefore offers a natural explanation for the presence of crystalline silicates in comets and the significant portions of crystalline silicates observed in accretion disks around young stellar objects.


1 Introduction

In the standard one-zone model of protoplanetary accretion disks (Pringle pri81 (); Lin & Papaloizou lin85 ()) the vertically averaged flow field of the disk is characterized by a radial inward drift of the disk matter within the inner, chemically active zone of the disk. Only in the icy region far away from the proto-sun the disk’s flow field is directed outward. However, Urpin (urp84 ()) found in his analytical work large-scale meridional flow patterns to exist. This meridional flow field is characterized by an outward directed drift of the disk matter close to the disk midplane whereas the flow is directed inward in higher layers. The result of Urpin (urp84 ()) has been confirmed by several authors by different analytical, semi-analytical and numerical methods (Siemiginowska sie88 (); Kley & Lin kle92 (); Różyczka et al. roz94 (); Kluźniak & Kita klu00 (); Regev & Gitelman reg02 (); Tscharnuter & Gail Tsa07 ()). Thus a meridional flow field seems to be the universal type of flow pattern existing in protoplanetary accretion disks.

With respect to radial mixing processes in protoplanetary disks the structure of the flow field is of utmost importance. This is a consequence of the fact that the average radial transport of matter in protoplanetary disks by advection occurs on a similar timescale as the transport by turbulent diffusion, namely on the viscous timescale. Hence, besides a realistic description of the turbulent diffusion, an exact knowledge of the large-scale flow field in disks is essential for calculating the transport of matter within disks. Here we consider the effect of meridional flows. Another type of hydrodynamic mixing associated with gravitational instabilities is discussed by Boss (bos04 (); bos07 (); bos08 ()) and found to be very efficient in that case.

Previous model calculations of protoplanetary disks which include the calculation of radial mixing of species have considered only the vertically averaged one-zone velocity field as the flow field of the disk (Stevenson & Lunine ste88 (); Cyr et al. cyr98 (); Drouart dro99 (); Bockelée-Morvan et al. boc02 () as well as the series of papers of the ITA group: Gail gai01 (); Wehrstedt & Gail weh02 (); Gail gai02 (); Gail gai04 (); Wehrstedt & Gail weh03 (), henceforth called Papers I – V). The only exception are the -dimensional model calculation of Keller & Gail (kel04 (), henceforth called Paper VI) and Tscharnuter & Gail (Tsa07 ()) where for the first time the meridional velocity field is used for the computation of the radial mixing of species in protoplanetary disks. The results show that the outward transport of species is much more efficient in the meridional flow field than in the one-zone flow field.

In the present work we extend the work of Paper VI. This is done by improving the time-dependent one-zone models of Papers II and V by calculating the disk structure in the +-dimensional approximation (e.g. Lin & Papaloizou lin85 ()) simultaneously with the -dimensional mixing of species in the disk. We therefore accept the small deviations of the +-dimensional disk structure from the disk structure of the exact -dimensional hydrodynamic calculations of Paper VI in order to save computing time. In contrast, however, the present +-dimensional disk model is coupled with a sophisticated chemical model which includes the calculation of equilibrium condensation of the most abundant solids, annealing of silicates, and combustion of solid carbon. Furthermore, the present disk model includes a detailed opacity calculation considering the Rosseland and Planck opacity means of the most important absorbers in the disk.

With this powerful tool we investigate radial mixing processes in protoplanetary disks. In particular the radial mixing of species in the solar nebula is of great interest to explain the composition of the most pristine solar system bodies, i.e., the comets. From observations it is long known that some fraction of the dust in comets is crystalline (e.g. Swamy et al. Swa88 (); Hanner et al. han97 ()) and crystalline silicate grains have been detected in accretion disks around young stellar objects (e.g. Meeus et al. Meu01 (); Bouwman et al. Bou01 (); van Boekel et al. vBo04 (); vBo05 (); Keller et al. Kel05 ()). Both these observations are interpreted as resulting from mixing material from the central parts of the disk to the outer regions though also other processes have been invoked for that (see, e.g., Alexander et al. Ale07 (); Wooden et al. woo07 () for a discussion).

The paper is organized as follows: Section 2 presents the treatment of the -dimensional transport of tracers in the disk within the disk model. In Sects. 3 and 4 the method of calculating the radial and vertical disk structure is described. Section 5 addresses the calculation of the velocity field of the disk, in particular the meridional flow field. Section 6 deals with the numerical treatment of the disk model. In Sect. 7 we present the results and finally make our conclusions in Sect. 8.

2 Transport of tracers

2.1 The transport-reaction equation

The transport of tracers with small concentration embedded in a carrier medium is given by the transport-reaction equation (Hirschfelder et al. hir64 (); cf. Paper VI)


where denotes the particle density of the carrier medium, the velocity vector of tracer component , the binary diffusion coefficient of tracer relative to the carrier, the rate term of gains and losses by chemical reactions of tracer , and


the concentration of tracer . denotes the particle density of tracer . The second term on the l.h.s. of Eq. (1) is the advection term. The first term on the r.h.s. the diffusion term.

By using the continuity equation


the transport-reaction equation (1) can be changed to


For simplicity


is assumed, i.e., all tracer velocities equal the velocity of the carrier medium. For the present model calculations this is an acceptable assumption since the considered tracers are micron-sized dust grains which are carried along with the flow in most parts of the disk. More precisely, for micron-sized particles Eq. (5) only holds in the densest and chemically active disk regions whereas in the less dense outskirts of the disk the motions of carrier gas and tracers may decouple, i.e., .

The binary diffusion coefficient is assumed to be given by


where is the kinematic viscosity of the disk matter and the Schmidt number of tracer . The Schmidt number is set in the present calculations to


for all tracers. The approximation for the Schmidt number (7) is based on the same approximations as assumption (5). A Schmidt number equal to unity means that the tracer moves as the carrier gas does. In less dense parts of the disk this approximation fails, and .

Since regions of low density do not contribute much to the total tracer transport, we apply (5) and (7) for all tracers in the present model. The value of can, however, not be fixed with any precision because of the unclear physics of the viscous transport in accretion disks. For some discussions on the value of see, e.g., Johansen et al. joh05 (); joh06 (); Turner et al. tur06 (); Pavlyuchenkov & Dullemond pav07 ().

For calculating the tracer transport in two dimensions we assume axial symmetry and use polar coordinates (). Here, denotes the polar radius and the polar angle measured from the midplane of the disk. With this choice of coordinates the transport-reaction equation (4) takes the form


where and are the velocities in radial and polar direction, respectively. The index denotes a kind of tracer occurring in a special modification (see following section).

2.2 The set of tracers

The tracers which are considered in the model calculations are:

  1. silicate dust grains (forsterite and enstatite) of different degrees of crystallisation and

  2. carbon dust grains of different sizes.

The tracers experience the following reactions:

  1. Silicate dust grains (forsterite and enstatite) start to anneal at temperatures above , i.e., the degree of crystallisation increases until it reaches unity.

  2. Carbon dust grains become decomposed by reactions with molecules at the grain’s surfaces. This process critically depends on the density of molecules in the gas phase and starts to operate at temperatures of under conditions encountered in protoplanetary disks.

The rate terms for the basic reactions are given in Paper II.

The set of transport-reaction equations (8) determines the concentration of crystalline forsterite, , and crystalline enstatite, , with different degrees of crystallisation , and the concentration of solid carbon grains with different grain sizes . The are used to calculate the average degrees of crystallisation of forsterite and enstatite (cf. Paper II)


and the degree of condensation of C into solid carbon


denotes the volume of a carbon atom within solid carbon and the (solar) abundance of C. is the number of sampling points for the discretised degrees of crystallisation , , and the number of sampling points for the size spectrum of solid carbon grains , . The set of sampling points for and is chosen as in Paper V (, ).

For a more detailed discussion of the processes of silicate annealing and carbon combustion the reader is referred to the other papers of this series, in particular Paper I and II.

2.3 Solution of the set of transport-reaction equations

The set of transport-reaction equations (8) is discretised in first order in time and in second order in space with respect to the diffusion terms. The advection terms are treated by a standard upwind method (see Paper II).

The transport-reaction equations (8) are solved by an ADI (Alternating Direction Implicit) method (Press et al. pre92 ()). By ADI each transport-reaction equation is split into two equations for separate directions in space: one equation contains the terms in -direction and the other equation the terms in -direction. The rate term is split into two half-steps and equally distributed on both directions. The resulting matrices of the separate equations have tri-diagonal structure and thus can be inverted numerically fast and easily (Press et al. pre92 ()). Inverting the matrix of the unsplitted transport-reaction equation (8), which has a band-tri-diagonal shape, would be numerically much more time-consuming and therefore less efficient than the ADI method.

To ensure numerical stability the sequence of solution of the splitted equations is permuted between two successive time steps (therefore Alternating Directions Implicit method). The transport-reaction equations are solved fully implicit.

3 Radial disk structure

The radial structure of the disk is calculated in the one-zone approximation by using the -prescription of viscosity and by accounting for the vertical self-gravity of the disk. The resulting set of equations for the radial disk structure is given by (cf. Paper V):

1.) Time evolution of the surface density :


where is the radial distance from the mass-centre.

2.) Keplerian angular velocity:


where is the Keplerian velocity, the gravitational constant and the stellar mass.

3.) -viscosity:


where denotes the viscosity parameter which is chosen to be in the model calculations (cf. Paper V).

4.) Isothermal sound speed:


where is the Boltzmann constant, the temperature in the midplane of the disk, the mean molecular weight and the proton mass, respectively.

5.) Pressure scale height by accounting for the disk’s vertical self-gravitation:


6.) Mean vertically averaged mass density:


7.) Mean molecular weight:


The calculation of the partial pressure of species in chemical equilibrium is described in Paper I.

8.) Rosseland and Planck means of the mass extinction coefficient:


For the calculation of the opacity the main dust absorbers are considered in the disk model. The abbreviations used here denote: ’for‘ (crystalline) forsterite, ’ens‘ (crystalline) enstatite, ’sil,am‘ amorphous silicate, ’car‘ solid carbon, ’iro‘ solid iron, ’cor‘ corundum and ’ice‘ water ice, ’mol‘ molecules. denotes the degree of condensation of the key element of condensate ( for the silicates, for solid carbon, for solid iron, for corundum and for water ice). The method of calculating , and is described in Paper I and II, respectively, and of , and in Paper V. and are calculated by Eqs. (9) and (10), respectively, and by Eq. (11). Analytical fit formulae for the Rosseland and Planck mean of the opacity of the individual species are given in Wehrstedt (weh03b ()).

9.) Rosseland and Planck mean vertical optical depth at the midplane:


10.) Viscous dissipation rate:


11.) Effective temperature of the disk surface:


where is the Stefan-Boltzmann constant and the temperature of the ambient molecular cloud.

12.) Temperature at the midplane:


13.) Radial drift velocity:


14.) Tracer transport in the one-zone approximation:


15.) Rate of mass accretion:


16.) Variation of the stellar mass by accretion of matter onto the protostar:


Note that is negative.

The set of equations for the radial disk structure (12) – (30) is solved with standard methods up to an accuracy of for the radial model structure.

4 Vertical disk structure

Since we wish to investigate the radial and vertical transport of tracers in protoplanetary disks as well as its feedback on the disk structure, it is necessary to resolve the vertical structure of the disk. For this purpose we choose the +-dimensional approximation for the calculation of the disk structure (e.g. Lin & Papaloizou lin85 ()). In this approximation the vertical stratification is calculated separately for each grid point of the radial one-zone model (see Sect. 3) by assuming a hydrostatic structure in the vertical direction. Therefore we accept some deviation of the results of the +-dimensional model calculation from that of the exact -dimensional hydrodynamic calculation, but avoid the numerical complexity of the latter.

4.1 Vertical disk equations

In the following the set of equations for the vertical disk structure in the +-dimensional approximation at a certain radius is given.

We define a -dependent column density:


where is the density at the height above the midplane at the given radius . Note that for with this definition. The value of is taken from the one-zone model. The differential form of Eq. (31) is


The pressure stratification is determined by


where denotes the pressure. The second term on the r.h.s. of Eq. (33) is the vertical acceleration of a self-gravitating infinite slab (Paczyński pac78 ()). The (vertical) self-gravity of the disk is taken into account since in Paper V it has been found that the self-gravity significantly modifies the disk structure for disk masses above . Since in the present model calculations we choose an initial disk mass of , the self-gravity of the disk has to be considered.

The generation of heat in the disk is assumed to be caused exclusively by viscous friction. Thus the vertical gradient of the energy flux is given by (e.g. Lin & Papaloizou lin85 ())


Note that holds according to Eq. (24). Further note that irradiation of the disk by the star is neglected in Eq. (34).

The solution of Eq. (34) requires the specification of the vertical dependency of the viscosity . However, the mechanism of turbulence in disks as the main source of viscosity is a matter of much debate so far (see e.g. Richard & Davis ric04 (); Johansen et al. joh05 (); joh06 (); Turner et al. tur06 (); Pavlyuchenkov & Dullemond pav07 () and references therein). Due to the lack of a precise knowledge of the vertical variation of , the viscosity at point (,) is set equal to the value which follows from the one-zone approximation at (Eq. (14)), i.e., the viscosity is set vertically constant.

The heat is assumed to be transported solely by radiation. Heat transport by convection is neglected since it is inefficient at low temperatures. In this case, the vertical gradient of the temperature is given by (e.g. Lin & Papaloizou lin85 ())


The Rosseland mean opacity at each vertical location is calculated as in the radial model (Eqs. (19) – (22)).

The set of vertical disk equations is closed by the equation of state for an ideal gas


4.2 Solution of the set of vertical disk equations

The set of equations for the vertical disk structure (32) – (36) constitutes a two-point boundary value problem. We define the top of the atmosphere at the height where the vertical optical depth is close to zero. According to the Eddington-Barbier approximation (e.g. Mihalas mih78 ()) the temperature at is determined by


The column density at has a very small value. In the model calculations we assume for a fixed value of


The energy flux at follows from Eq. (34) to be


Finally the pressure at has to be determined. By assuming an isothermal stratification for small optical depths Eq. (33) can be integrated. Its solution is



denotes the error function and the mean molecular weight at which is set to .111 The temperature at is sufficiently low that the matter is in molecular form at any radial position at the top of the atmosphere, and hence is close to for the standard cosmic element mixture.

Starting with an initial guess of the height , the set of vertical disk equations (32) – (36) is solved by a Runge-Kutta method of th order (Press et al. pre92 ()) from down to the midplane of the disk at . The Runge-Kutta solver is equipped with a step size control which keeps the number of vertical grid points small. In general, the column density at found by this procedure is not equal to the midplane value . Then is iterated by a shooting method until within an accuracy of the vertical model of .

5 Velocity Field

Finally, the velocity field in the -dimensional disk model has to be specified. We will show in this paper that the flow field in the disk affects the radial mixing of tracers in the disk, exceeding the contribution of turbulent diffusion. Therefore the actual flow field in the protoplanetary disk has to be determined.

5.1 Velocity field of the one-zone approximation

The velocity field of the one-zone approximation is given by Eq. (27). The main feature of this velocity field is that there exists a characteristic radial position where the sign of the radial velocity changes:

  • Inward of , is negative, i.e., the flow is directed inward.

  • Outward of , is positive, i.e., the flow is directed outward.

In model calculations of the solar nebula is typically located at a few or a few tens of and slowly changes with time (e.g. Ruden & Lin rud86 (); Paper II). This is far outside of the chemical active zone of the disk which is located inside of about . Hence in one-zone models of the solar nebula, outward radial mixing of tracers from the chemical active zone to the outer solar system is only possible by turbulent diffusion (Paper II; Paper V).

5.2 Meridional velocity field

The real velocity field in disks is more complex than those of the one-zone approximation. Urpin (urp84 ()) considered higher order terms in the flow equations of the disk and found a characteristic flow pattern:

  • close to the disk midplane, the flow is directed outward whereas

  • at elevated heights above the midplane the flow is directed inward.

This type of meridional flow field has been confirmed to exist by -dimensional hydrodynamic calculations by Kley & Lin (kle92 ()), Różyczka et al. (roz94 ()), Tscharnuter & Gail (Tsa07 ()). Regev & Gitelman (reg02 ()) conclude that the meridional flow field is a dynamical phenomenon that is of universal occurrence in accretion disks.222In the -dimensional extension the flow field of accretion disks should additionally show large vortices induced by the baroclinic instability (Klahr & Bodenheimer kla03 (); Klahr kla04 ()).

In the case of disks with -viscosity, Takeuchi & Lin (tak02 ()) derived analytic expressions for the radial and azimuthal velocity and , respectively. In their work was calculated in first order of the small perturbation parameter . is calculated in second order of (which is the lowest non-vanishing order). Keller (kel03 ()) and Paper VI obtained the same result as Takeuchi & Lin (tak02 ()) and additionally obtained the vertical velocity in st order of .

In the following the meridional velocity field for disks with -viscosity is derived in the lowest non-vanishing order of . According to Takeuchi & Lin (tak02 ()) and Keller (kel03 ()) the velocity field in the stationary limit is given by


As mentioned above, the lowest order, no-vanishing deviation of from Keplerian rotation is of second order in . as well as are of first order in . To simplify Eqs. (41) – (43), the density distribution of the stationary +-dimensional disk model is used, i.e. (e.g. Pringle pri81 ())333 In Eq. (44) the self-gravity of the disk is neglected. It would be an interesting task to include the disk’s self-gravity in the calculation of the meridional velocity field.




is the midplane density and


the scale height. denotes the radius of the protostar. Finally a radial power law for the midplane temperature is assumed,


where is the radial temperature exponent444 In principle, this assumption is not needed. The quantity in all equations for the velocity can simply be interpreted as the local value of d / d . The impact of the vertical temperature structure on the meridional flow field would be also of interest. However, this is not considered in the present work. . For these density and temperature profiles, the velocity field in -disks is as follows:


with the factor

and coefficients

To get a first impression of the meridional flow field see Fig. 4b.

As in case of disks with -viscosity (Takeuchi & Lin tak02 (); Paper VI) for the meridional velocity field of -disks (Eqs. (48) – (50)) there exists a characteristic height where the radial velocity changes its sign:


provided the quantity under the square root is positive, which either requires (if ) or . The latter is not realised in accretion disks models calculated with reasonable dust opacities.

If , which is satisfied in most parts of the disk (see Sect. 7.3.1 and Fig. 4(b)), is negative above , i.e., the flow is directed inward, whereas is positive below , i.e., the flow is directed outward. Otherwise, in the case , is negative at any height, and the disk matter (the mixture of gas and small particles) moves inward at all heights.

The factor affects the flow pattern only in the very inner parts of the disk (). As a consequence the location of , and therefore the flow pattern of the disk, essentially depends on the radial temperature exponent .

In the model calculations of Paper VI, has been chosen as a fixed parameter and set constant throughout the disk. In the present work is calculated locally, i.e., at each grid point from the radial temperature profile by a -point interpolation. The interpolated value of is used to calculate the -dimensional velocity field of the disk (,) given by Eqs. (49), (50).

A vector transformation is implemented in the code for transforming the -dimensional velocity field (,) into polar coordinates (,) since polar coordinates are used to calculate the transport of tracers (cf. Eq. (8)).

Note int Eq. (48) that the angular velocity deviates from the Keplerian velocity . In particular it varies with height . However, the value of from Eq. (48) is not required in the present -dimensional model.

6 Numerical treatment

6.1 Initial and boundary conditions

The boundary conditions are chosen as follows:

  • At the outer boundary () we set (no-torque condition).

  • At the inner boundary () we choose the quasi-stationary boundary condition which has been introduced in Paper V. The quasi-stationary boundary condition removes the problem of the unphysical behaviour of in the inner zone of the disk which appears in the case of the no-torque condition and ensures a smooth and physically more realistic radial distribution of the surface density (as well as other quantities) in the inner disk zone.

  • The concentrations at the outer boundary are chosen such that the silicates are amorphous () and all carbon that is not bound in is condensed into solid carbon ().

  • The at the inner boundary are chosen such that the silicates are completely crystalline () and solid carbon is completely oxidised ().

  • At the upper boundary of the -dimensional polar grid () the dust is assumed to have interstellar properties as is assumed for the outer boundary. Hence, the at the upper boundary are chosen such that the silicates are amorphous () and a fraction of of the total carbon is assumed to be bound in solid carbon grains, the remaining fraction being in .

  • For the upper boundary as well as the inner and outer boundary homogeneous von Neumann boundary conditions are chosen for the tracres, i.e., we set at the upper boundary and at both lateral boundaries. This choice prevents any flux of matter across the boundaries.

  • At the midplane of the disk symmetry conditions are applied.

The initial conditions for the one-zone model are chosen as in Paper V for the model with solar abundance, i.e., the initial radial distributions of and are calculated from a stationary model ( and ) in which the silicates are assumed to be unaltered amorphous ISM grains.

The at the position (,) of the -dimensional polar grid initially are set equal to the midplane value , i.e., the initially are set constant in -direction.

Parameters (, , , )Initial model (, )Set of radial disk Eqs. (13) – (27)NR iteration of (,) up to Solution of Eq. (12) for Solution of the set of Eqs. (28) for the New , and Fixpoint iteration of up to Set of vertical disk Eqs. (32) – (36)Iterating up to an accuracy of for Interpolation from the (,) grid onto the gridSolution of the set of Eqs. (8) for the Interpolation from the (’,) grid onto the gridNew , and Fixpoint iteration of up to New stellar mass from Eqs. (29) and (30)Next time step
Figure 1: Flow chart of the code. The arrows to the right of the diagram denote the iteration loops and the time loop, respectively. The inner two loops in the upper half refer to the one-zone model and the inner two loops in the lower half to the -dimensional model. denotes the accuracy of the iteration. ‘NR’ is an abbreviation for Newton-Raphson. The set of radial disk Eqs. (13) – (27) is solved by a coupled Newton-Raphson method for the midplane temperature and the mean molecular weight . For details see text.

6.2 2-dimensional grids

In the present model the disk structure is calculated in cylindrical coordinates (,) (+-dimensional approximation) whereas the tracer transport is calculated in polar coordinates (,). For a self-consistent computation of the tracer transport and the disk structure an interpolation of the relevant physical quantities between both grids is applied. This is done by a D bilinear interpolation (Press et al. pre92 ()).

The radial grid extends from to and consists of logarithmically equidistant grid points ( grid points per decade). The vertical grid size depends on the radial position and is variable in time since we adopt a step size control for in the calculation of the vertical disk structure. In the chemically active zone of the disk the number of vertical grid points usually is larger than (maximum grid points) whereas it drops below in the cool outer disk regions. The total number of grid points of the cylindrically symmetric grid adds up to about .

For the polar grid an apex angle of has been chosen. This choice ensures that the polar grid is embedded in the cylindrically symmetric grid at any time of the model calculations. The polar angle is discretised in equidistant grid points between the upper boundary and the midplane . Note that the polar grid defines the ’domain of transport‘, i.e., the domain where the transport of tracers by advection and diffusion occurs during the model calculations. Beyond the polar grid no transport takes place (, ). We did test calculations with polar grid points as well as with an apex angle of and found no significant deviations from the standard model with and .

6.3 Numerical Solution

The flow chart of the code of the present disk model is shown in Fig. 1. The model calculations are performed as follows:

First the radial disk structure in the one-zone approximation is calculated with standard methods (inner two loops in the upper half of Fig. 1; see Sect. 3 and cf. Papers II and V). The equation for , Eq. (12), is solved fully implicit.

Secondly, the vertical disk structure is computed (innermost loop in the lower half of Fig. 1; see Sect. 4). In this way we obtain the disk structure in cylindrical coordinates (,).

Thirdly, the -dimensional transport of tracers is calculated. This is done by

  1. Interpolating relevant quantities (, , , , , ) from the cylindrically symmetric grid (,) onto the polar grid ,

  2. solving the set of D tracer equations (8) in polar coordinates (see Sect. 2),

  3. calculating the degrees of crystallisation of the silicates, and , as well as the degree of condensation of in solid carbon, , from the (see Eqs. (9) – (11)), and

  4. interpolating , and from the polar grid onto the cylindrically symmetric grid (,).

To calculate the tracer transport self-consistently with the disk structure the temperature is iterated globally by a fixpoint iteration up to an accuracy of (outer loop in the lower half of Fig. 1). The relative low accuracy of the global iteration is attributed to the errors which are introduced by the interpolation procedure. We did test calculations without performing the global iteration and found only small deviations from the standard model with which are not relevant for the final results.

Finally, the stellar mass is updated (see Eqs. (29) and (30)) before the next time step is executed (outermost loop in Fig. 1).

6.4 Model parameters

Initial stellar mass
Stellar effective temperature
Stellar luminosity
Stellar radius
Inner disk radius
=0.096 AU
Outer disk radius
Molecular cloud temperature
Initial disk mass
Disk angular momentum
Apex angle of the polar grid
Viscosity parameter
Number of polar
 grid points =
Number of cylindrical
 grid points
Table 1: Parameters used for the calculation of the disk models.

The model parameters for the disk models calculated in the present work are shown in Table 1. The stellar and disk parameters are chosen as in Paper V and are typical for solar like T Tauri stars and their surrounding accretion disks (see references given in Paper II). With respect to the choice of the value of the viscosity parameter the reader also is referred to Paper V and the discussion therein.

The model calculations are initiated with a small time step of which is increased slowly by an implemented time step control. The time step is limited to to ensure numerical stability. We quit the model calculations at . The model with the above standard parameters requires a CPU time of about two month on a P4 XEON computer.

7 Results

7.1 Models 1DC and 1DP

Before we present the results of the D model calculations we first compare two model calculations in the one-zone approximation which correspond to two different underlying disk geometries.

In the first model, henceforth called model 1DC, the tracer transport is calculated in cylindrical coordinates (,) by assuming vanishing concentration gradients in the vertical direction, i.e. , . The transport-reaction equation in the one-zone approximation then is given by Eq. (28), i.e.,


Model 1DC essentially equals the standard model with solar element mixture in Paper V.

In contrast to model 1DC, the tracer transport in the second model, henceforth called model 1DP, is calculated in polar coordinates (,) within a one-zone model. For this purpose we assume the concentration gradients in polar direction to be negligible, i.e., . This assumption is almost identical to the assumption made for Eq. (52). It is justified by the fact that the disk model is based on the idea of a flat and flaring disk () which has a small apex angle. As well it is justified by the results of the D model calculations of the present work which show only slight concentration gradients in the vertical (respectively polar) direction since the vertical (polar) diffusion erases almost any vertical concentration gradient.

With these basic settings the transport-reaction equation in D polar coordinates, Eq. (8), transforms into


Equation (53) determines the tracer transport in the one-zone model 1DP.

Both Eqs. (52) and (53) determine the tracer concentrations in the midplane of the disk where but the geometry of the disk in both models is different. In model 1DC the disk is a flat slab whereas in model 1DP the disk resembles an outward flaring slab. For this reason in both models the diffusion terms in the transport-reaction equations, Eqs. (52) and (53), differ from each other.

Figure 2: Tracer transport in one-zone models with different geometries at times (0), (1) and (2). The ’cylindrical‘ model 1DC is shown with dashed lines and the ’polar‘ model 1DP with solid lines. (a) Degree of crystallisation of forsterite . (b) Fraction of condensed in solid carbon .

The radial velocity profile in both models is assumed to be that of the one-zone model which is given by Eq. (27).

Note that the radial density profiles that enters in the diffusion terms of Eqs. (53) and (53) in both models is identically to the density profile of a flaring disk (cf. Eq. (17)). However, to compare the tracer transport in models 1DC and 1DP one actually has to apply a density profile for model 1DC with constant scale height to account for the cylindrical disk geometry in this model. In contrast, the polar disk geometry in model 1DP represents to be a good approximation of the real flaring disk geometry. We realised this inconsistency after the termination of the model calculations. However, we refrain from recalculating model 1DC with a density profile with constant scale height since we intend to compare the model calculations of the present paper with those of the previous papers of this series which computes the tracer transport as model 1DC, i.e., with cylindrical disk geometry and ’flaring‘ .

The dependence of the results for the tracer transport on the different assumed disk geometries is shown in Fig. 2.

In Fig. 2a the degree of crystallisation of forsterite at (0), (1) and (2) is plotted for the models 1DC (dashed lines) and 1DP (solid lines). It clearly can be seen that the ’polar‘ model 1DP leads to a substantially less efficient radial mixing than the ’cylindrical‘ model 1DC. In the region around after of the disk evolution, in model 1DP is about one order of magnitude lower than in model 1DC. This can be explained by the disk geometry. In model 1DP the disk flares. As a result in model 1DP tracers are diluted to a greater extent as they are mixed outward as compared to model 1DC where the geometry is similar to a tube with parallel walls. However, polar coordinates obviously resemble more closely the real disk geometry than cylindrical coordinates do. Thus, in former one-zone models of the solar nebula that treat the tracer transport in cylindrical geometry (Paper II; Paper V; Bockelée-Morvan et al. boc02 ()) radial mixing has been overestimated. As a consequence the hypothesis of radial mixing that explains many interesting properties of primordial solar system bodies hardly can be maintained for previous one-zone models. For example the observed fraction of crystalline silicates of more than of the bulk silicate in many comets (e.g. Hanner et al. han94 (); Crovisier et al. cro97 (); Hanner et al. han97 (); Yanamandra-Fisher & Hanner yan99 (); Wooden et al. woo05 (); woo07 ()) hardly can be explained by one-zone models with a realistic polar disk geometry.

A way out of this dilemma is provided by -dimensional models with a realistic description of the flow field in the disk (see Sect. 7.3).

Figure 3: +-dimensional grid of the standard model at in the inner part of the disk. Each cross marks a grid point. The dashed line shows the radial profile of the pressure scale height (Eq. (16)) and the solid line the upper boundary of the D polar grid. The big dots close to the upper boundary of the grid mark the location of the photosphere (). The arc-like structures mark the location of the vapourisation zones of strong absorbers (for details see text).

Carbon dust in model 1DP is mixed radially outward more efficiently than in model 1DC (Fig. 2b). This can be explained by the dilution of the products of carbon combustion in the outer parts of the disk (e.g. , ; Gail gai02 ()) which is more pronounced in the polar model 1DP than in the cylindrical model 1DC. This again is due to the geometric effect described above for the radial mixing of the crystalline silicate grains.

7.2 1+1-dimensional grid

In the next step we turn from the one-zone disk model to the +-dimensional disk model with -dimensional transport of tracers.

To give an impression of the +-dimensional cylindrically grid it is plotted in Fig. 3 at time of the standard model 2DM for the innermost of the disk and . As been described in Sect. 4.2, the solver for the disk vertical structure is equipped with an automatic step size control keeping the number of vertical grid points small. In particular the step size is small in such regions where physical quantities show large gradients. Such regions clearly can be seen in Fig. 3.

On the one hand the grid point density gets high close to the upper boundary of the vertical grid as a consequence of the large pressure and density gradients.

On the other hand the grid point density is increased in the narrow zones of dust vapourisation/destruction since temperature and opacity are strongly variable there. As a result, arc-like structures of enlarged grid point density are formed which stand out from the ambient grid with low resolution. An example is the arc-like structure that culminates at height and extends out to . This results from the transformation of enstatite to forsterite. Above and to the right of the arc enstatite and forsterite are both stable whereas beneath the arc enstatite is unstable and forsterite is the only stable silicate. The enstatite-to-forsterite transformation occurs at about . The two arcs culminating at and display the vapourisation fronts of forsterite () and solid iron (), respectively. The vapourisation of corundum which is the most refractive dust species in the present model calculations occurs at and (midplane values at ).

For comparison, in Fig. 3 also the radial run of the pressure scale height (Eq. (16); dashed line), the upper boundary of the D polar grid (solid line) and the location where the atmosphere gets vertically optically thin (; big dots) is shown. The D polar grid for tracer transport (not shown) is well embedded in the +-dimensional cylindrically symmetric grid. Moreover, the D polar grid is well located within the vertically optically thick region of the disk.

Note again that the D polar grid defines the area of turbulent and advective transport. Beyond these area no diffusive or advective transport is considered in the present model calculations.

7.3 2-dimensional model calculations

In the final step we present the results of the 2-dimensional model calculations. For the purpose of investigating the influence of the disk‘s flow field on the tracer transport two model calculations with different flow patterns are performed:

  • In the first model, henceforth called model 1DE, the flow field of the one-zone model (Eq. (27)) is extended to the disk regions off the midplane, i.e., at each polar grid point the polar radial velocity is chosen equal to the one-zone radial velocity , , and the polar velocity everywhere is set to zero, . Although this velocity field is -dimensional we call it in the following the ‘extended one-zone’ velocity field (1DE).

  • In the second model, henceforth called model 2DM, the meridional velocity field is used which is calculated as described in Sect. 5.2.

7.3.1 Flow fields in models 1DE and 2DM

To compare the extended one-zone velocity field and the meridional velocity field the flow fields of both models 1DE and 2DM are plotted in Fig. 4 at time in the region between and .

The one-zone velocity field (Fig. 4a) is directed inward for most parts of the disk as it was mentioned in Sect. 5.1. The velocity increases inward as the viscous torque increases with decreasing . The absolute value of the velocity at , e.g., is whereas at the inner boundary it is . The feature with large grid point density at about marks the ice condensation front.

Figure 4: -dimensional flow field in the disk at between and . The solid line marks the upper boundary of the polar grid (). (a) Model 1DE. (b) Model 2DM. The arrows show the direction of the flow; their length is proportional to the velocity. For details see text.

The meridional flow field (Fig. 4b) shows a distinctly more complex structure than the one-zone velocity field. The meridional flow pattern shows large eddies which extend from both sides of the midplane of the disk to the disk’s surfaces. Particularly the structure of the flow field in the region of the ice condensation front is of interest as the radial temperature distribution shows a complex behaviour in that region. In detail the midplane radial velocity shows the following radial dependence in model 2DM at :

  1. For , is positive since for the radial temperature exponent in this region holds (cf. Sect. 5.2).

  2. In the range of , is negative. There the opacity is a decreasing function of that results in a radial temperature profile with .

  3. In the range , again is positive as the condensation front of ice leads to a temperature plateau, i.e., .

  4. In the range , for the same reasons as in (ii) is negative.

  5. Finally, for , is positive since the temperature adjusts to the ambient cloud temperature (), i.e., .

Due to the accretion process the regions (i) – (v) slowly move radially inward. Additionally, the diffusive spreading of the disk causes a flattening of the radial temperature profile which narrows the ranges (ii) and (iv) with negative with time. Despite of this, the flow patterns (i) – (v) around the ice front persists for more than of the disk evolution. The chemical active zone where in particular annealing of silicates and carbon combustion occurs, is always located within range (i).

The polar velocity is more than one order of magnitude lower than close to the midplane, and therefore plays no important role with regard to mixing processes in the disk.

Figure 5: Tracer transport in the -dimensional models 1DE (dashed lines) and 2DM (solid lines) at times (0), (1) and (2). Also shown is the result of the polar one-zone model 1DP (dashed dotted lines). (a) Degree of crystallisation of forsterite . (b) Fraction of condensed in solid carbon .

7.3.2 Radial mixing: models 1DP and 1DE

We now study the impact of the meridional flow field on the efficiency of radial mixing processes in the disk. For this purpose in Fig. 5 the midplane radial distributions (a) of the degree of crystallisation of forsterite and (b) of the fraction of condensed into solid carbon are displayed for the models 1DE (dashed lines) and 2DM (solid lines) with the same representation as in Fig. 2. For a comparison the result of the polar one-zone model 1DP is plotted in Fig. 5 with dotted lines.

Figure 6: Time evolution of the degrees of crystallisation of forsterite, , (filled symbols) and enstatite, , (open symbols) in the models 1DE (squares) and 2DM (circles) at the position of the ice front.

We first compare the models 1DP and 1DE. In both models the tracer transport is calculated in polar geometry with the same one-zone velocity field. Hence these models show the influence on the tracer transport if extending it from D to D. The main difference of models 1DP and 1DE can be attributed to their different temperature structures. In the present work the models in the +-dimensional approximation generally lead to lower midplane temperatures for the chemically active parts of the disk than models in the one-zone approximation. This is an explicit result of the calculation of the vertical disk structure in the +-models, particularly of calculating the radiative transfer in vertical direction. The midplane temperature difference between the models 1DP and 1DE can amount to more than in the innermost parts of the disk. As a consequence chemical reactions and dust processing in model 1DE occur at a smaller distance than in model 1DP, i.e., the annealing front of silicates (Fig. 5a) as well as the carbon combustion front (Fig. 5b) in model 1DE are shifted inward as compared to model 1DP.

Model 1DE Model 2DM

Table 2: Degrees of crystallisation of forsterite and enstatite, and , respectively, in the model calculations 1DE and 2DM at some selected times and radial positions .
Model 1DE Model 2DM