First simulation results of Titan ’ s atmosphere dynamics with a global 3-D non-hydrostatic circulation model

We present the first results of a 3-D General Circulation Model of Titan’s atmosphere which differs from traditional models in that the hydrostatic equation is not used and all three components of the neutral gas velocity are obtained from the numerical solution of the Navier–Stokes equation. The current version of our GCM is, however, a simplified version, as it uses a predescribed temperature field in the model region thereby avoiding the complex simulation of radiative transfer based on the energy equation. We present the first simulation results and compare them to the results of existing GCMs and direct wind observations. The wind speeds obtained from our GCM correspond well with data obtained during the Huygens probe descent through Titan’s atmosphere. We interpret the most unexpected feature of these data which consist of the presence of a non-monotonicity of the altitude profile of the zonal wind speed between 60 and 75 km.


Introduction
Saturn's largest satellite Titan is the only moon in the solar system that possesses a considerable atmosphere with a surface pressure about 1.5 times that of the Earth (e.g.Taylor and Coustenis, 1998).Titan's synchronous rotation period of barely 16 days implies that, like Venus, Saturn's largest Correspondence to: H. K. Biernat (helfried.biernat@oeaw.ac.at) satellite belongs to the group of slowly rotating bodies with a dense atmosphere.Under such conditions, the circulation of the atmosphere is not close to the geostrophic balance, as on Earth or Mars, where pressure gradients are balanced by Coriolis forces, but by the centrifugal forces of an almost latitudinal circulation.In Titan's case the atmosphere is in a cyclostrophic balance (Flasar, 1998;Grieger, 2004).Titan's atmosphere has been subject to several studies with Global Circulation Models (GCMs), providing results that differ considerable depending on the models applied (e.g.Müller-Wodarg at el., 2000;Tokano et al., 1999;Hourdin et al., 1995).This may be due to the fact that the temporal evolution of the temperature, which is used to derive the horizontal and vertical wind dynamics in a GCM, is mainly governed by advective and diffusive heat transports and the radiative heating term.The latter one is strongly influenced by both the aerosol absorption properties and their spatial distribution, which, however (up to now), are not very well constrained.In turn, the aerosol distribution may also be influenced by the atmospheric circulation (e.g.Cabane et al., 1992) which implies a coupled system.
Nowadays, the widely used GCMs for the Earth's atmosphere, as well as the GCMs developed for other planetary atmospheres, are mainly based on the hydrostatic equation (e.g.Roble at el., 1988;Wilson and Hamilton, 1996;Forget at el., 1999;Müller-Wodarg at el., 2000).In our study we introduce a new mathematical model, which provides the three-dimensional global distributions of the zonal, meridional, and vertical wind components without any restrictions on the vertical transport of the neutral gas.
In the current version of our model we do not implement a radiative transfer model to derive the atmospheric temperature profile and its diurnal changes, but assume a predefined Published by Copernicus GmbH on behalf of the European Geosciences Union.and stationary temperature profile.This was adapted from the currently available Titan atmosphere model based on Voyager 2 ultraviolet, infrared, and radio occultation measurements (Yelle et al., 1997).This assumption clearly reduces the complexity of the GCM, as it avoids the necessity to model the planet's aerosol properties and distribution, and the coupling of haze formation with atmospheric dynamics, as presented by Rannou et al. (2002).It should be noted that the approach of a GCM with a predescribed relaxation temperature field was successfully demonstrated for Titan's tropospheric circulation by Grieger (2004).
In Sect. 2 we outline the mathematical formulation of our GCM, it's computational implementation, and the assumptions used in the simulations.In Sect. 3 we briefly review the current knowledge of Titan's atmospheric dynamics that was derived from the currently available ground and spacebased observations.In Sect. 4 we describe our simulation results and compare them to the currently available observational data.In Sect. 5 we summarize the current status of the model, identify shortcomings and physical inconsistencies arising from simplifying assumptions, and provide an outlook on the possibilities for future applications of our model.
We note that the successful application of our model to the Earth's atmosphere was presented in Mingalev andMingalev (2001, 2005).In particular, the model was applied to investigate mechanisms responsible for the formation of the circumpolar vortices at levels of the stratosphere and mesosphere in the Northern and Southern Hemispheres.The vortices, obtained from our model, were found to be consistent with the global circulation of the atmosphere, obtained from observations.

Coordinate system, model region, and grid spacing
We use a Titan centered rotating Eulerian spherical coordinate system with the radial (vertical) coordinate r being the distance of the current position to the planet's center, and the azimuthal and poloidal coordinates designated as the longitude ϕ and the latitude β, respectively.Our model covers a 250-km thick spherical layer above the planet's surface.The longitudinal and latitudinal grids are uniform, with a step size of ϕ= β=1 o .The radial grid above 72 km has a constant step size of 2 km.Below 72 km the radial grid is nonuniform, with a step size that slowly decreases from 2 km down to 0.255 km at Titan's surface.

Fundamental equations
The version of our model presented enables us to calculate the 3-D global distributions of the gas density ρ, as well as the vertical, zonal, and meridional wind components, v r , v ϕ and v β , respectively.Positive directions for the wind components are defined as upward, eastward, and northward.The three-dimensional global distributions of the gas density, as well as the vertical, zonal, and meridional wind components, are obtained from the numerical solution of a system of equations comprising the continuity equation the vertical wind component the zonal wind component (3) and finally, the meridional wind component Equations ( 2)-( 4) are the components of the dynamical equation (Batchelor, 1970).Furthermore, we take into account the ideal gas law, given by In Eqs.
(2)-( 5) p denotes the gas pressure, T the temperature (with a defined distribution), R the universal gas constant, and M the mean molecular weight, which was adapted from the Y97 recommended Titan atmosphere model (Yelle et al., 1997).=4.6×10 −6 s −1 is the angular velocity of Titan's rotation, g 0 =1.35 m/s 2 is the gravitational acceleration on the planet's surface, r 0 =2575 km is Titan's radius, and rr , ββ , ϕϕ , rϕ , rβ , ϕβ are the physical components of the viscous stress tensor.The last tensor is defined by and takes into account the turbulent friction.Here I is the identity tensor and Tr( D) is the trace of the tensor D, which is defined as where v is the wind velocity vector, ∇⊗v is the tensor of gradient v , (∇⊗v) T is the transposed tensor of gradient v , η is the symmetric tensor of the viscosity coefficients, which can be expressed in the given model coordinate system by where µ r =µ+ρ K r , µ ϕ =µ+ρ K ϕ , and µ β =µ+ρ K β .Here µ=1×10 −5 kg m −1 s −1 is the molecular viscosity coefficient taken from Müller-Wodarg et al. (2000) and K r , K ϕ , and K β are the turbulent viscosity coefficients in the vertical, zonal, and meridional directions, respectively.In our model these coefficients are given by the Richardson formula: where C=0.37 m 2/3 s −1 (Richardson, 1926) and L r , L ϕ , L β are the vertical, zonal, and meridional scales of turbulence.
We define these scales as follows: where h, ϕ, and β are, respectively, the radial, longitude, and latitude grid step sizes.Our simulation results suggest an insignificant sensitivity to variations of the turbulence scales up to 50% in either direction.
Substituting Eqs. ( 8) and (7) into Eq.( 6), one obtains the components of the viscous stress tensor: Because of Titan's negligible magnetic environment (Rishbeth et al., 2000), no ion-neutral coupling was taken into account in Eqs. ( 2 The numerical solution of Eqs. ( 1)-( 5) is obtained from a special version of an implicit splitting difference scheme, which was adopted for each equation of the system.To be able to control the relaxation of the numerical solution to the stationary solution, each of Eqs. ( 1)-( 4) is implemented in the form where A i are the terms of the equation.Furthermore, we calculate for each of Eqs. ( 1)-( 4) the nonstationarity parameter ε , (at each point of grid), which is defined as It is assumed that the numerical solution sufficiently approximates the stationary solution, if the maximum value of ε on the grid for all Eqs.( 1)-( 4) does not exceed 0.03.

Test of the model
To check our model we consider the class of exact solutions of the system (1)-( 5), which is given if the temperature depends only on the potential U (i.e. the sum of the centrifugal and gravitational force).In this case the wind is zero and the gas density depends only on U where U o is the value of U at the point where the prescribed value of the gas density ρ is known.Using this temperature dependent distribution, a zero initial wind, and an initial density distribution (which does not satisfy Eq. 9) the model provides a solution after the relaxation process, which approximates the exact stationary solution with a satisfactory accuracy.The numerical solution obtained possesses little fluctuation in the vicinity of the immobility.The amplitudes of the fluctuations were less than 10 −4 m/s.This fact implies that we can calculate the stationary circulation of Titan's atmosphere, corresponding to the given temperature field by using our model.

Titan's atmosphere dynamics from observations
An observational basis for understanding Titan's atmosphere dynamics and meteorology was virtually nonexistent prior to the Voyager 1 flyby on 12 November 1980.In particular, valuable information was obtained from infrared (Hanel et al., 1981) and radio science (Tyler et al., 1981) observations.Current reviews on Titan's dynamic meteorology is provided (among other authors) by Bird et al. (2005) and Flasar (1998).Strong zonal winds of ∼100 m/s are implied by the latitudinal temperature gradient deduced from the Voyager 1 infrared observations (Flasar et al., 1981;Flasar et al., 1992).Assuming a hydrostatic, gradient-balanced flow, the zonal wind velocity can be related to the latitudinal temperature gradient by the "thermal wind equation".If one neglects the winds near ground level, this equation may be solved for the zonal wind height profile, provided one has knowledge of the temperature with height and latitude.Under these conditions, however, the thermal wind equation does not constrain the wind direction, i.e. whether the winds are prograde or retrograde.Voyager 1 provided the following equator-to-pole temperature gradients observed at three levels: T 16 K at the 0.5 mbar pressure level at ∼230 km from observations in the 1304 cm −1 channel (Flasar et al., 1981;Coustenis, 1990); T 1 K at p=100 mbar (∼40 km) from observations at 200 cm −1 ; and T 2 K at the surface from the thermal channel at 530 cm −1 .Toon et al. (1988) have cautioned, however, that these observations may include a significant contribution from stratospheric aerosols.If the IR-observations are stretched to the maximum possible equator-to-pole gradients, corresponding to double the bestfit values (Lunine et al., 1991;Flasar et al., 1997), a zonal wind height profile can be derived with zonal velocities a factor of the order of √ 2 larger than the best estimate.Figure 1 (a reprint of Fig. 1 of Bird et al., 2002) shows a summary of what the observations and previous GCM simulations imply about Titan's zonal winds.The zonal wind height profile is plotted at latitude of 19 • N. The latitude-adjusted profile derived from the Voyager 1 observations by Flasar et al. (1997) is shown by the short-dashed curve in Fig. 1.The solid line is the corresponding "maximum" envelope for the given latitude.The dash-dotted curve is a mean of the northern spring equinox models of Hourdin et al. (1995).Figure 1 also shows linear fits to the results of a Titan dynamical process study using a terrestrial GCM, incorporating various idealizations of the stratification and drag likely to be relevant to the planet's specific circulation regime (Del Genio et al., 1993; longdashed lines T1-T4).
Infrared heterodyne observations of Titan's ethane emission in the 12-µm band, which originates near the 1-mbar level (∼200-km altitude), have been used to estimate Titan's zonal wind speeds (Kostiuk et al., 2001).The results provide strong evidence that Titan's zonal winds are moving prograde (with a 94% probability derived from three separate observation opportunities at the NASA/IRTF at Mauna Kea).The wind speeds derived from simple models for the hemispheric mean measurements are stated as 210±150 m s −1 .
The traditional cloud-tracking technique to derive winds was not very successful from Voyager 1 images.The Titan encounter by the Cassini spacecraft on 2 July 2004 provided the first opportunity during this mission to measure cloudtracked winds on Titan by using the Visual and Infrared Mapping Spectrometer (VIMS) (Momary et al., 2004;Griffith 2002)).The curve given by the short-dashed line is the integral of the thermal wind equation, using equator-to-pole contrasts at various heights inferred from the Voyager 1 infrared observation (Flasar et al., 1997).The solid curve is a "maximum envelope" (Lunine et al., 1991), assuming upper bounds on the Voyager temperature gradients.The dash-dotted curve is a prediction from a specifically developed Titan circulation model (Hordin et al., 1995).The long-dashed curves labeled "T1-T4" are linear fits to four atmospheric simulations for a planet rotating at the Titan period (16 days) performed by Del Genio et al. ( 1993) over a more limited range of altitudes: T1: model with global absorbing cloud at altitudes from 20-40 km; T2: same atmospheric simulation, but without an absorbing cloud; T3: T1 without stratospheric drag; T4: T1 with greatly weakened surface drag.et al., 2004).Spectral imagery revealed that cloud coverage of Titan was sparse, covering less than 1.5% of the observed sunlit surface, nevertheless, several clouds were followed during the encounter.The most prominent cloud was located near the South Pole (roughly circular with a diameter of 600±110 km) and was tracked over 11 images spanning a 13-h period.Momary et al. (2004) derived a prograde mean wind speed of 0.5±3.3m/s.
The results from the Huygens descent were published by Tomasko et al. (2005).The results from the Huygens Doppler Wind Experiment were published by Bird et al. (2005).The altitude profile of Titan's zonal wind was derived from Green Bank Telescope (GBT) in West Virginia and Parkes Radio Telescope in Australia throughout the course of the Huygens probe descent through the atmosphere of Titan.The Huygens probe entered Titan's atmosphere at the point with spatial coordinates of 10.33±0.17• S latitude, 196.08±0.25 • W longitude, and 154.8±11.2km altitude, began to descend, and landed on Titan after approximately 150 min.Over the duration of the descent, the Huygens probe drifted eastward following the horizontal wind.The estimated longitude of the Huygens landing site on Titan is approximately 192.33±0.31• W. Figure 2 (a reprint of Fig. 2 of Bird et al., 2005) shows the variation of the zonal wind with altitude and pressure level derived from GBT and Parkes ground-based Doppler data.The most striking feature of the measured profile is the region of strong reversed shear between the 65 and 75-km altitude, where the zonal wind decreases to a minimum of 4 m/s, which then reverts to a strong prograde shear above 75 km.This wind shear layer was unexpected and is unexplained up to now.

Discussion of results and comparison with observations
We make an attempt to simulate numerically the steady-state global distributions of the horizontal and vertical wind in Titan's atmosphere for conditions characteristic for the Huygens probe descent.It is convenient to present the simulation results in a specific spherical coordinate system.This coordinate system is obtained from the Titan planetographic coordinates by the following transformation.The centers of both coordinate systems are assumed to be identical and situated at the center of Titan.The principal axis of the spherical coordinate system is turned so that it remains in the plane comprising both the vector of Titan's angular velocity, , and radius vector from the center of Titan to the center of the Sun, r o .Moreover, it is supposed that the principal axis of the spherical coordinate system becomes perpendicular to the vector r o and is directed toward the Northern Hemisphere of Titan.It is obvious that the sub-solar point on Titan is located at the equator of the new spherical coordinate system, that is, at the line where the latitude is equal to zero, i.e. β=0.
It is assumed that the longitude of the sub-solar point on Titan did not change during the transformation of the coordinate system and has a value of 156.85 • W. In new coordinate system, the point of atmospheric entry of the Huygens probe is situated at approximately 7.3 • N latitude and 164 • E longitude.(Flasar et al., 1997).The estimated uncertainties in the zonal wind speed, based on an adaptation of an error analysis for the Huygens-Cassini link, are of the order of 80 cm/s at high altitude and drop roughly in proportion to the absolute speed to 15 cm/s just above the surface.These uncertainties are primarily systematic errors associated with the Huygens trajectory at entry.The estimated statistical (measurement) error is always smaller, the standard deviation being of the order of <5 cm/s towards the end of descent.With the possible exception of the region above 100 km, where the wind fluctuations are greatest, the zonal flow is found to be generally weaker than those of the model.The wind shear layer in the height range between 60 and beyond 100 km was unexpected and is at present unexplained.
We take into account the fact that the Titan's angular velocity is rather slow (16 times less than that of the Earth).A period of revolution of Saturn around the Sun is approximately 30 times as long as the Earth's year.Consequently, it may be expected that the temporal evolution of the global temperature and wind distributions in Titan's atmosphere are much less than those in the Earth's atmosphere.Therefore, in model calculations, we assume that the global temperature distribution is stationary and conditioned by solar illumination, with the spatial structure of the atmospheric temperature being invariable relative to the sub-solar point on Titan.Hence, the global temperature field is symmetric relative the equator of the specific spherical coordinate system utilized in the model calculations.It can be noticed that the equator of the specific spherical coordinate system coincides with the equator of the planetografic coordinates at the equinox sharp.Thus, in the model calculations, the spatial structure of atmospheric temperature is analogous to the global temperature distribution characteristic for the equinox on Titan.

Utilized temperature distribution in Titan's atmosphere for equinox
We approximate Titan's troposphere, stratosphere, and mesosphere experimental temperature field for equinox conditions by the distribution T (h, ϕ, β), which is given by the formula In Eq. ( 10), h denotes the altitude, T r (h) denotes the vertical temperature profile (taken from the Y97 recommended model and shown in Fig. 3.), ϕ o denotes the longitude of noon meridian (taken as 156.85 • W in our calculations), T ad (h) describes a vertical temperature profile, which is used to adjust the temperature field for latitudinal variations and makes it consistent with observational data, and T z (h) describes a vertical temperature profile, which is used to adjust the temperature field for longitudinal variations.
In the 0-45-km altitude region the profile T ad (h) where X= [(h 1 −h)/ h 1 ] 2 and h 1 =45 km and in the 45-180-km altitude region is given in [K] by formula In the 180-250-km altitude region From Fig. 3 and from Eq. ( 10) one can see that the pole-toequator (more precisely the intersection of the noon meridian and the equator) temperature differences are 16 K above 180 km, and 2 K on Titan's surface, and therefore are consistent with the Voyager 1 equator-to-pole gradients, as mentioned in Sect.3. Equation ( 10) furthermore implies that the constructed temperature distribution is characterized by a difference between the noon meridian -equator intersection and the midnight meridian -equator intersection temperature of 0 K on Titan's surface and of only 1 K at an altitude of 250 km.

Boundary conditions and initial conditions for the model calculations
As was noted in Sect.3, different experimental data manifest that, in Titan's atmosphere at heights near 200-km altitude,  Kostiuk et al. (2001).On Titan's surface, all components of the wind speed vector are assumed to be zero.These boundary conditions ensure the conservation of mass for the entire model region.
Our assumption is zero initial meridional and vertical wind components in the model region and a horizontally homogeneous initial gas density, as given by the Y97 recommended model.Furthermore, we adapted the Y97 recommended model is (spatially homogeneous) chemical atmosphere composition.The initial zonal wind component in the model region is given in [m/s] by formula v ϕ (h, β)=260 (h/ h o ) cos β , where h o =250 km and h is the altitude in km.

Discussion of the simulation results
For the simulation of the stationary circulation of Titan's atmosphere we used the global temperature field given by Eq. ( 10) and the initial conditions as described in the previous section.
Figures 4a-c  Figures 6a-c show the corresponding plots for the vertical wind component for the respective altitudes of 40, 120, and 200 km.
The zonal, meridional and vertical wind components distributions for the longitude of 164 • E are depicted in Figs.7a-c, respectively.
It is doubtless that the calculated circulation of Titan's atmosphere is conditioned by the horizontal nonuniformity of the temperature and strong superrotation at the upper boundary.This circulation possesses the following features.The flow of Titan's atmosphere is symmetric relative to the equator.Throughout the wide height range attached to the upper boundary, the horizontal wind velocity is primarily zonal and eastward.In the lower atmosphere, a significant meridional component of the wind velocity arises at middle latitudes.Near the equator, from 5 km down to the surface, the horizontal wind velocity is predominantly directed to the equator and can have a west zonal component.This feature is in accordance with data obtained during the Huygens probe descent through the atmosphere of Titan.This west zonal component of the wind velocity can achieve a magnitude of a few m/s.
The results of the simulation indicate that the horizontal and vertical components of the wind velocity are changeable  functions not only of latitude but also of longitude.The dependence of the calculated functions on longitude is more pronounced at higher altitudes.
It can be seen from Fig. 7a, in particular, that a large-scale irregularity of zonal wind distribution exists at the vertical plane lain along the meridian, containing the point of atmospheric entry of the Huygens probe, over the latitude range from 20 • S to 20 • N at altitudes below 100 km.
The altitude profiles of the zonal wind velocity at the points near the equator ought to possess the following features.Above 30 km the calculated value of zonal wind velocity rises up to an altitude of 60 km, where it achieves a maximum, then decreases up to 75 km, where it achieves a minimum, and then rises monotonically as far as the upper boundary.Figure 8 shows the zonal wind height profile derived from our GCM (solid line) at a latitude of 7.3 • N and a  longitude of 164 • E which correspond to the point of atmospheric entry coordinates of Huygens probe.It can be noticed that the altitude profile, derived from our GCM (solid line in Fig. 8), differs from the profiles obtained throughout the course of the previous GCM simulation and shown in Fig. 1.
Figure 7b shows that, at the longitude of 164 • E in the 45 • S-45 • N latitude region, the meridional wind component in the altitude ranges from 0 to 12 km and from 35 to 65 km is directed from the poles towards the equator and in the altitude ranges from 15 to 30 km and from 70 to 250 km is directed from the equator towards the poles.

Comparison with observations
The 200 m/s zonal wind speeds in the 200-250-km altitude region correspond very well with the observed upper atmospheric prograde wind measurements from Kostiuk et al. (2001) in the 0.1-0.7-mbarpressure range.
The upward vertical mass transport in the pole regions at levels of the middle and upper troposphere in our simulation results fit well together with the reported observations of cloud formations in Titan's upper atmosphere's south pole region (see Sect. 3).This flow could produce cloud formations near the poles due to condensation of CH 4 .Comparing the GCM results in Fig. 8 to the results from the Huygens Doppler Wind Experiment, one can see that the solid line in Fig. 8 is qualitatively close to the Huygens profile.Both profiles possess the following features.The regions of strong reversed shear between 60 and 75-km altitude are present, where the zonal wind speed decreases to a minimum, which then reverts to strong prograde shear above 75 km.Over the height range from 45 to 125 km, the difference between the observed and calculated values of zonal wind speed does not exceed 20 m/s.Thus, model calculations have permitted us to simulate the most unexpected feature of the zonal wind speed profile, measured during the Huygens probe descent (see Lebreton et al., 2005) through the atmosphere of Titan, which consists of the presence of the non-monotonicity of this profile, no doubt caused by a specific global circulation of the Titan's atmosphere, characterized by the presence of largescale, three-dimensional vortices.The meridional motion, as well as vertical mass transport, taken into account by the global non-hydrostatic circulation model utilized, plays a significant role in the formation of these vortices.The calculated three-dimensional circulation of Titan's atmosphere is affected by the strong superrotation at the upper boundary and is analogous to the flow between two concentric spheres revolved with different angular velocities.In Titan's atmosphere, in addition, this flow is influenced by the presence of the horizontal nonuniformity of the temperature.As a consequence, near the equator, the narrow band arises, where the large-scale irregularity of zonal wind distribution exists, which produces the non-monotonicity of the altitude profile of the zonal wind velocity between 60 and 75 km.

Present status of the model and simulation results
The main characteristics of the current version of the GCM can be summarized as follows: -The use of full Navier Stokes equations, allowing for non-hydrostatic flows.
-An imposed temperature profile, which was derived from observations and is considered as stationary.
-A strong super-rotation at the upper boundary of the model region, which corresponds with experimental results.
-A very high horizontal and vertical resolution.
The main features of the simulation results of Titan's atmospheric dynamics can be summarized as follows: -In the polar regions, meridional and zonal winds are weak.
-In the 45 • S-45 • N latitude region, the horizontal mass transport changes direction in the altitudes range of 0-100 km.The meridional wind component in the altitude ranges from 0 to 12 km and from 35 to 65 km is directed from the poles towards the equator and in the altitude ranges from 15 to 30 km and from 70 to 250 km is directed from the equator towards the poles.Also, the meridional wind attains a speed of about 20 m/s at altitudes of 100-150 km, about 8 m/s at altitudes of 50-60 km and 20-35 km, and about 12 m/s in the altitude range of 0.5-4 km.
-Values of the vertical winds lie in the range of -1.5 to 0.7 m/s and may be higher than in reality (see further discussion in the subsequent section).
Above 100 km, the zonal wind velocity is eastward and rises with the altitude to 150-200 m/s near the equator.Below 100 km, the zonal wind velocity is eastward, too, and, over the latitude range from 20 • S to 20 • N, possess the nonmonotonicity of the altitude profile between 60 and 75 km.This feature corresponds very well with data obtained during the Huygens probe descent through the atmosphere of Titan.The non-monotonicity of the altitude profile of the zonal wind speed between 60 and 75 km is undoubtedly caused by the meridional, as well as vertical mass transport, creating the specific large-scale, 3-D vortices over the wide height range near the equator, where from 5 km down to the surface, the horizontal wind velocity has the west zonal component which is consistent with the measurements carried out with the help of the Huygens probe.
It may be recalled that the cloud coverage of Titan was observed experimentally in the polar regions (see Sect. 3).The simulation results indicate that the vertical wind velocity is upward below 60 km in the polar regions.It may be expected that the simulated vertical mass transport could promote cloud formation due to condensation of CH 4 near the poles.

Problems and future work
The main shortcoming of the current model version is the use of an imposed temperature profile and the assumption of stationarity.This is especially important, as one has to keep in mind that the adiabatic heating by radiation drives the circulation through pressure gradient forces and the temperature distribution is the result of the whole system.The lack of energy equations produces a problem in which the vertical winds obtained in our current simulations may be stronger than expected.A direct constraint on the mean vertical wind arises from the adiabatic heating/cooling associated with vertical subsidence/ascent.The estimation of the adiabatic heating rate can be made very easily from the model results.For instance, vertical winds of the order of 3 m/s associated with vertical gradients of potential temperature of typically 1 K/km (g/C p =1.3 K/km on Titan, so that an isothermal vertical temperature profile corresponds to a vertical gradient of potential temperature of 1.3 K/km) will produce heating rates of the order of several K per hour.This heating or cooling can only be balanced by radiation, on average, in the real atmosphere.The problem is that the expected radiative heating by solar and thermal infrared radiation is by order of magnitude smaller than this, except maybe at the top of the model (Flasar et al., 1981).As a consequence, there is no term in the real Titan atmosphere which can balance the huge cooling associated with the vertical ascent predicted by the model at the equator.If the model would include an energy conservation equation rather than imposed temperature, the equatorial temperature would rapidly decrease at the equator and increase at the pole with such a meridional circulation.
This work clearly shows that one of the main necessary updates of our model is the implementation of the energy equation by a proper radiation transfer model which incorporates microphysical processes as haze condensation and sublimation.
The peculiar dynamics of the atmospheres of slowly rotating planetary bodies, characterized by a coherent zonal motion and super-rotation in low latitudes, has been subject to theoretical studies with GCMs for a long time.It is now generally accepted that the main driver for super-rotation is the thermal tide, the energy of which is converted into regular zonal motion due to weak Coriolis forces and hence the low amplitude of Rossby waves.As thermal tides are substantially non-balanced modes of atmospheric motion, we expect that our non-hydrostatic model will more adequately describe the related phenomena than the traditional GCMs.

Fig. 1 .
Fig. 1.Models of zonal wind height profiles on Titan at a latitude of 19 • N (reprinted Fig. 1 from Bird et al. (2002)).The curve given by the short-dashed line is the integral of the thermal wind equation, using equator-to-pole contrasts at various heights inferred from the Voyager 1 infrared observation(Flasar et al., 1997).The solid curve is a "maximum envelope"(Lunine et al., 1991), assuming upper bounds on the Voyager temperature gradients.The dash-dotted curve is a prediction from a specifically developed Titan circulation model(Hordin et  al., 1995).The long-dashed curves labeled "T1-T4" are linear fits to four atmospheric simulations for a planet rotating at the Titan period (16 days) performed by Del Genio et al. (1993) over a more limited range of altitudes: T1: model with global absorbing cloud at altitudes from 20-40 km; T2: same atmospheric simulation, but without an absorbing cloud; T3: T1 without stratospheric drag; T4: T1 with greatly weakened surface drag.

Fig. 2 .
Fig. 2. Titan zonal wind height profile (reprinted Fig. 2 from Bird et al. (2005)).The zonal wind derived from GBT and Parkes observations is compared with the prograde Titan engineering wind model and envelopes based on Voyager temperature data(Flasar et al., 1997).The estimated uncertainties in the zonal wind speed, based on an adaptation of an error analysis for the Huygens-Cassini link, are of the order of 80 cm/s at high altitude and drop roughly in proportion to the absolute speed to 15 cm/s just above the surface.These uncertainties are primarily systematic errors associated with the Huygens trajectory at entry.The estimated statistical (measurement) error is always smaller, the standard deviation being of the order of <5 cm/s towards the end of descent.With the possible exception of the region above 100 km, where the wind fluctuations are greatest, the zonal flow is found to be generally weaker than those of the model.The wind shear layer in the height range between 60 and beyond 100 km was unexpected and is at present unexplained.

Fig. 3 Fig. 3 .
Fig. 3 and 5a-c show the distribution of the horizontal wind component for the respective altitudes of 1, 40, 80, 120, 160, and 200 km.The color code provides the wind speeds in [m/s].
Figure 7b also shows that the meridional wind attains a speed of about 20 m/s at the altitudes of 100-150 km, about 8 m/s at the altitudes of 50-60 km and 20-35 km, and about 12 m/s in the altitude range of 0.5-4 km.

Figures 6
Figures6 and 7cshow that the vertical wind attains a speed of about -1.5 m/s and a speed of about 0.7 m/s.Figures 6a and 7c also show that, in the polar regions, vertical wind is directed upward in the altitude range of 0-55 km.Figures4, 5, 7a, and 7b show that in the polar regions meridional and zonal wind components are small.Figures4a and 7ashow that there is no superrotation at altitudes of 0-4 km.Figure4bshows that zonal wind component at the altitude of 40 km is directed eastward and attains values of more than 40 m/s in the 35 • S-20 • S and 20 • N-35 • N latitude re-

Fig. 8 .
Fig. 8.The dashed line shows the zonal wind height profile derived from the Huygens Doppler Wind Experiment.The solid line depicts the zonal wind height profile derived from our GCM at a latitude of 7.3 • N and a longitude of 164 • E which corresponds to the point of atmospheric entry of the Huygens probe.