PAPER II

Introduction

In this paper, a detailed theoretical study is made to determine the exact terms in the tangent-linear model (TL) and its adjoint (AD) which are responsible of the mutual influence existing between prognostic and diagnostic model variables in a four-dimensional variational analysis. These terms are then quantified in experimental investigations aiming at the determination of the synoptic ingredients responsible for some abnormal features like the over/under-estimation of specific humidity increments.

The first section is dedicated to the nonlinear, tangent-linear and adjoint theoretical formulations. The second interprets some experimental investigations performed with complete and partial sets of observations and tries to link the interpretations to the expressions found for the above-mentioned terms.

1. Eulerian tangent-linear formulation of the primitive model equations

To understand the multi-variate effects in a 4d-var assimilation system, which uses the model dynamics and physics as information sources (among others) to find the best linear unbiased estimate (BLUE), the hydrostatic primitive equations of ARPEGE/IFS are described in a humid adiabatic context. Then a linearized form of these equations at the vicinity of a basic nonlinear model state is presented; this will underline the different terms that govern the behaviour and magnitude of the gradients of the different control variables. This will show, in particular, the terms influencing the specific humidity adjoint calculations.

The framework is assumed as Cartesian. This greatly simplifies the form of horizontal spatial operators (as e.g. horizontal derivatives). But, this does not change at all the essence of all the features discussed in this paper, compared to the way they appear in the real (non Cartesian) model. The vertical coordinate is the same as in the numerical model; it is the hybrid pressure coordinate h ( h Î[0,1]), which is defined implicitly from two (arbitrary) functions A and B by :

where P is the hydrostatic pressure and P _s the surface hydrostatic pressure. The following boundary conditions are imposed for A and B :

1.1 The nonlinear system equations

Two versions of the equations are used for the dynamics in ARPEGE/IFS : the Eulerian version and the semi-Lagrangian one. All the experiments described below are achieved in a 4d-var setup using Eulerian formulation both for trajectory and TL/AD integrations. That's why all the equations are expressed in Eulerian form.

The Cartesian humid adiabatic system for Eulerian 3d rotating atmosphere directly reflects the corresponding nonlinear equations in the numerical model :

where T is the temperature, the vector W represents the earth rotation, the horizontal wind vector, f the geopotential, q the specific humidity, R the air perfect-gas constant and Cp the specific heat at constant pressure. "Ñ " and "Ñ×" stand respectively for the horizontal gradient and divergence operators.

We have also the following diagnostic relationships :

The first one can be obtained by vertical integration of the hydrostatic equation; the second one is a form of the continuity equation, while the two last ones are obtained by vertical integration of the surface pressure equation.

The pressure-force term in the momentum equation " RTÑln(P)+Ñf " includes a contribution due to the horizontal pressure gradient and another one given by the geopotential gradient. This term is very important in the discussions below. It depends on almost all the prognostic variables (q via the air perfect-gas constant R =(1-q )R_d + qR_v where R_v and R _d are respectively the water vapour and the dry air constants, T, P and P_s ), it implies a vertical integration and an horizontal pressure gradient. It constitutes a dynamic term allowing multi-variate mutual influences for these fields. It will be studied in details when formulating its linearized expression.

Another important term, with a similar impact, is the conversion term which appears in the thermodynamic equation "RT  /C_p ·w  /P". R and Cp  =  (1-q)Cp_d  + qCp _v (where Cp_d and Cp_v  are the specific-heat constants for dry air and water vapour respectively ) depend both on specific humidity; w is a diagnostic function of wind components and pressure. It appears, then, that this term is also responsible for mutual influence between variables during a 4d-var assimilation process.

The other terms which are mainly the horizontal and vertical advections and Coriolis force play also the same role, but in a less important way.

1.2 Linearized form of the primitive Eulerian hydrostatic equations.

There can be several different forms for the linearized versions of the unique original hydrostatic primitive adiabatic system. This comes from the fact that the basic state can be chosen in several ways. In the context of 4d-var analysis, the linearization is done in the vicinity of a basic state called "trajectory" obtained by integrating the nonlinear model in parallel with its linear integration. This linearization integration is then called the tangent-linear one.

When linearizing model equations, one assumes that all terms implying the products of perturbations are very small, and then could be neglected. This is correct only if the time integration is short and the model trajectory doesn't cover a rapid baroclinic atmospheric situation.

The prognostic variables used in the TL/AD versions of the forward nonlinear model are the perturbations u', T', q ' and P_s' of momentum, temperature, specific humidity and surface pressure. The linearized equations will be written in a form implying only these variables; the diagnostic parameters such as f, w and y will be linearized first, and expressed using the prognostic perturbations.

The RT term which depends on both temperature and specific humidity will be considered as one variable renamed Z. We have : Z'=(R_v  - R _d)Tq'+RT'

The nonlinear hydrostatic equation takes the expression ¶_hf =Z· ¶_h ln(P), with horizontal gradient :

The linearized form of this equation is then :

The only boundary-condition required for integrating vertically this equation is the horizontal gradient of the surface-geopotential perturbation which is assumed to be zero. This leads to :

On the other hand, if we assume : ¶ _hP =m, then : m'=¶ _hP'. The linearized expression of is a function of ' and P'. After this change of variable, it reads :

The linearized momentum equation can then read :

The last term includes contributions of both q' and T'. It is equal to the following expression :

This form shows a linear evolution of ' depending on ', P', T' and q' as :

The linearized thermodynamic equation reads :

The last term in this equation could be divided in two parts depending respectively on q' alone and T' alone as :

Thus the linear evolution of T' is expressed as :

The linearized equation for specific humidity is simpler as it consists only in linearizing the corresponding horizontal and vertical advections terms :

This may be written as :

Finally, the last equation to linearize is the evolution of the pressure field. The 3d-pressure perturbation may be deduced from the perturbation of surface pressure by : P'=B(h ) P'_s . The linearized form of the equation for pressure is then :

The following system summarises the four prognostic linearized equations :

By transposing the matrix appearing on the right-hand-side of this system, we obtain the adjoint form of the linearized equations.

When neglecting the temporal time-step (or forcing it to one second - 1.-), the specific humidity gradient is given by the following equation :

When we assume dry TL/AD integrations (and the same for the corresponding trajectory) by considering that R  = R _d (or R _v  =   R_d ), we eliminate the impact of C_u and we simplify C _T .

In the ARPEGE/IFS setups, Cp_v and Cp_d are set as :

Cp_v   =  4.   R_v and Cp _d  =  3.5  R_d .

Then,

Cp_v   =  8./7.   Cp_d and Cp   =  (7.+ q)/2.  R   =  (7.+ q)/ 7.   Cp _d ,

R/Cp  = 2./( 7.+q) »  2./7. and (Cp_v -Cp_d)/Cp  =  1/(7.+q) »  1./7. ,

so that :

RT /C _p·w  /P·[ (R_v  - R _d)/R  - ( Cp_v-Cp_d )/Cp] is reduced to : -2.Tw  /P  /(7.+q) ² »   -2./49.   Tw / P .

But, when we put R _v  =  R _d at the level of CVA1 (under key LDRYTL ), the setups are already done. In this case :

RT /C _p·w  /P·[ (R_v  - R _d)/R  - ( Cp_v-Cp_d )/Cp] is reduced to R_d Tw  /P·(4.R_v  -3.5 R_d)/[4. qR_v  + 3.5  (1-q )R_d] ² .

The effect of C_T is still present but, due to its sign, it affects the specific humidity increments in the opposite way.

In the real case, where R =(1-q)R _d + qR_v and Cp  =   (1-q) Cp_d  + qCp _v we have :

where q is considered to be dimensionless. That is, we have : (R_v  -  R _d)/R   £    (Cp_v -Cp_d )/Cp .

To have an idea about the trajectory parameters that contributes in the anomaly, we can compare the values of C_T before and after applying R _v  =   R_d . If we denote them C_T(bef) and C_T(aft) respectively, then at the beginning of the assimilation window we have :

This difference is large in an atmosphere characterized by high temperatures, strong vertical motions, low pressures. This context is typical of highly baroclinic atmosphere.

The adjoint model integration is a backward one having as initial state the last tangent-linear forecasted perturbation on the assimilation window, transported to the observation space thanks to the observation operators, then normalized by the observation errors and returned back to the physical space by the adjoint of the observation operators, as could be seen on the analytical expression of the cost-function gradient :

where B is the matrix of model error covariances, R_i the matrix of observation error covariances for time-slot "i", H' and M' the tangent-linear versions of the observation operator and the nonlinear model respectively (^T indicates the transposition of the matrixes), and d_i is the innovation vector.

The process of TL/AD integrations is iterated several times, but the trajectory is calculated only once. That means that the coefficients of the matrixes involved in the TL/AD models are computed once, whereas the perturbations dx are iterated several times as required by the minimizing process.

2. Experiments achieved to go further in the investigations of the specific humidity analysis anomaly for the situation of April 24th, 2001, at 18:00 UTC

2.1 Description of the experiments

The first experiments, described in Paper I, consisted in determining which observation types contributed the most to the anomaly. The aim was to identify minimum degrees of freedom in order to be able to interpret easily the difficult 4d-var assimilation processes involved in an experiment. This sensitivity study showed a great impact of AIREP and SYNOP types. Each of these two types is responsible for about 50% of the large negative specific humidity increments (Figs  2a-b and 5a).

More precisely a single SYNOP observation of surface geopotential at 18:00 UTC is found to be responsible for about 30% of the anomaly caused by all the SYNOPs. The corresponding station is an automatic one, identified as 08233, localized at [40.93°N,-1.30°W]. Two AIREPs, EU4593 and EU1456, are also concerned.

The experiments below will use this single SYNOP geopotential observation to diagnose the anomaly and to quantify the several terms involved in the computation of the specific-humidity increments. The strategy that will be followed consists in quantifying the magnitudes of each term (C _u, C_T , C_q), and then deduce in which context this problem could happen knowing the ingredients (in the trajectory or in the perturbations) that are responsible for the magnitude of the main term among these three.

This will lead to more visibility on the context causing this kind of anomaly. Thus, we could identify and handle efficiently the real cause of the problem.

2.2 Increments from a single observation : theoretical aspects

In case of a single observation, situated at time-slot i >  0, the incremental cost-function has the following form :

Here d_i    =  D and   R _i^-1  =  s₀ ^-2 are reduced to scalar values.

The analysis increment dx _a is obtained for Ñ J =  0. After some rearrangements, and if we suppose that the observation is situated at the gridpoint corresponding to the nth element of the analysis vector, then it is expressed as :

It is now clear that the analysis increments at the beginning of the assimilation window are obtained from the normalized analysis increments at time-slot i, transported by the adjoint model to time-slot 0, and then normalized by the model errors.

This preliminary theoretical study shows that if one wants to analyse the behaviour of TL/AD in 4d-var in the context of a single observation, this later must not be situated at the beginning of the assimilation window. In this case only observation operators would act in the minimization, and the problem would be a 3d one.

In the experiments described here, the observation is situated at the middle of the assimilation window and the increments are evaluated and visualized at the beginning of this window. Let's remember that the abnormal drying introduced by the operational analysis in this situation concerns all the time-slots between 15:00 and 21:00  UTC.

If we analyse the increments at the middle of the window, their expression is given by :

This case is similar to a 3d-var case, where the background cost-function (distance to the guess) implies the matrix MBM ^T. It is then clear, that in presence of rapidly evolving atmospheric phenomena the model errors may be greatly amplified by the tangent-linear model. Thus, assimilation will rely more on the observation ingredient. This will be the only information governing the assimilation.

The absence of specific-humidity observations combined with the above-mentioned fact could lead to abnormal analysis structures.

2.3 Contribution of SYNOPs and AIREPs to the specific-humidity increments.

The parameters from SYNOPs used in the upperair analysis are 10-m wind and geopotential. The other surface parameters are not used. The observed surface-geopotential gives information on temperature and pressure via the observation operators. So, introducing a SYNOP geopotential is similar to using an observation of temperature and pressure. The signal on temperature is then converted to a specific-humidity increment by the strong multi-variate property of 4d-var (Fig.  5).

The parameters used from AIREP reports are horizontal wind and temperature. In this situation, AIREPs were very dense on Western Europe along the assimilation window, as shown on Fig.  1.

2.4 Impact of AIREP reports on specific-humidity increments for April 24th, 2001, at 15:00

An experiment consisting in running 4d-var assimilation with only AIREP reports showed that this observation type contributed a lot in creating large specific-humidity increments (up to 15.  10^-4) along a band including the anomaly zone (Fig. 2a). AIREPs caused also another area of large specific-humidity increments (up to 13.  10^-4 ) North-East of the anomaly zone.

When combining TEMP and AIREP reports into one same experiment, the North-East part of the band of large specific-humidity increments disappears entirely. This fact is associated to the presence in this region of 3 TEMPs as shown in Fig.  3b. Let us remember that along the assimilation interval there is no TEMP reports in the entire anomaly area over Spain. TEMP reports are the only observation types informing on upperair humidity.

Figure 2 : 850 hPa specific-humidity increments at 15:00 UTC, obtained using : a) AIREPs only b) a complete set of observation types. The values are multiplied by 10⁴.

On the other hand, unexpectedly, the AIREP parameter associated to specific-humidity increments anomaly is not temperature, but wind (Figs 4a-b). As can be seen from the equations presented in part 1.2, the term acting here is C_u^· u ^ad .

Figure 4 : Contribution of AIREP data to the anomaly : a) temperature, b) winds . The values are multiplied by 10⁴ .

2.5 Impact of SYNOPs on specific-humidity increments for the situation of April 24th, 2001

An experiment consisting in using only SYNOPs in the cost function has been performed. As shown on Fig.  5a, the impact on specific-humidity increments is also large (about 15. 10^{- 4}) and corresponds to the anomaly zone. The parameter responsible for that is mainly the geopotential (by 99 %); 10-m wind has a far smaller impact (Fig.  5).

Figure 5 : Contribution of SYNOP information to the specific-humidity increments : a) all data, b) only geopotential data . The values are multiplied by 10⁴. Fields at 850 hPa are shown, as for the the other plots of increments.

Since 10-m wind has a neutral effect on humidity analysis, and surface geopotential acts as an observation of temperature, the adjoint term that is involved in the amplification of humidity increments here is mainly C_T  ^·T^ad .

These diagnostic experiments show that the anomaly is caused by two separate ingredients associated to AIREPs and SYNOPs. In the following, the impacts of single-observation experiments involving these two types are shown.

2.6 Impact of a single SYNOP geopotential observation on specific-humidity increments.

The following table gives the characteristics of the concerned SYNOP geopotential observation :

Figure 6 : Initial trajectory fields at 850 hPa for : a) temperature, b) vertical velocity, c) pressure, d) geopotential. These fields are the same in all experiments, they correspond to the guess fields at low resolution.

Figure 7 : Increments at 15:00 UTC caused by a single geopotential observation for : a) specific humidity, b) geopotential, c) temperature, d) vertical velocity, .

When looking at the vertical zonal and meridional slices of increments for the different prognostic variables, we realize that all the information comes from the surface (Fig.  8), indicating the vertical influence of a surface geopotential observation.

Figure 6a shows a temperature field with relatively large values in a thermal dorsal regime (values between 4°C and 14°C) over Spain, associated to large "RT" horizontal gradients (Fig.  9d), large geopotential gradients (Fig. 6d), relatively high vertical velocities (Fig.  6b) and low pressures (Fig.  6c).

As said in part 1.2, the trajectory quantity that contributes to the specific-humidity increments (when an observation informing on temperature is used) is mainly :

Figure 9c shows the corresponding field, demonstrating that when an atmospheric situation presents simultaneously the following 3 ingredients, specific humidity has a great probability to behave in an abnormal way over a region devoid of humidity observations :

1. High temperatures
2. Low pressures
3. Great vertical velocities

The representation of temperature, vertical velocity and pressure on Figs 6a-c show that these fields verify the above-mentioned conditions over Spain leading to a large difference in C_T(bef)   -  C_T(aft) »   0.17Tw   /P . This term is represented in Fig.   9c and shows effectively absolute minima over the entire domain. The term C_T · Grad(T) gives a minimum over Spain corresponding very well to the anomaly area.

Figure 8 : Vertical slice on a zonal axis of increments at 15:00 UTC for : a) specific humidity, b) temperature, c) geopotential, d) vertical velocity .

Figure 9 : Derived fields at 850 hPa and 15:00 UTC : a) C_T multiplied by 10⁵ , b) C_T ·Grad(T), c) C _T(bef)  - C _T(aft) »  0.17 Tw /P, multiplied by 10 ⁵, d) RT .

2.7 Impact of a single AIREP wind observation on specific-humidity increments.

The following table gives the characteristics of the concerned AIREP wind observation :

It is not possible to find one single AIREP wind contributing alone significantly to specific-humidity increments. All individual observations contribute in the same proportion. That is why the figures represented below (Fig. 10) show a little (but not negligible) effect of wind data from a single AIREP element. But if we use all the observations available in one AIREP report (as EU4593 or EU1456), the effect is then visible and corresponds to the summation of all the individual effects of each elementary observation (Fig.  11). The term involved here is obviously :

The trajectory quantities governing its magnitude are the temperature and its horizontal gradient, the horizontal and vertical gradients of pressure. If we define the operator n as :

then this term has the following form :

Figure 10 : (a) Specific-humidity increments at 850 hPa caused by one AIREP wind (the minimum is 6.10^{- 5}), (b)  corresponding temperature increments at 850 hPa, (c) vertical slice of specific-humidity increments at latitude 37° N, (d) geopotential increments in the same environment.

Figure 11 : Specific-humidity increments at 850 hPa caused by : a) EU1456 report and b) EU4593 report. (c) and (d) are the corresponding vertical slices.

3. Conclusion and discussion

One problem of the humidity analysis is the appearance, in a certain atmospheric context, of large analysis increments at the lowest model levels, over some subtropical land areas and also sometimes over temperate regions. These areas are characterized by high temperatures and dry conditions. Investigations done by Anderson et al., 1998, with ECMWF 3d-var, showed that these humidity increments were caused by geopotential observations rather than humidity data. Geopotential data can be fitted by changing both temperature and humidity. In the absence of any other data, the relative changes of humidity and temperature when fitting geopotential data are governed by the background-error standard deviations (of temperature and humidity). Anderson et al. concluded that there was a problem with the specification of humidity background-errors in hot and dry conditions. An easier but less correct solution would have been to disable the dependency on humidity of the geopotential observation operator.

This was a problem observed in 3d-var data assimilation. With 4d-var this problem is still present and more complicated to diagnose, because in 4d-var the observation operators include in addition the forecast model. There can also be a great influence on the humidity analysis from surface pressure and upperair geopotential or temperature data, through at least the virtual-temperature effects of the pressure force (including its geopotential part and its pressure part), and the conversion term of the thermodynamic equation : in addition to the relations involved in observation operators.

Rabier et al., 1998, showed that 4d-var analyses are moister than 3d-var analyses, by up to 0.3 g/kg when zonally averaged along the equator. The interpretation given to this fact is that the effective background-error standard deviations are larger than for 3d-var, taking into account the fact that the implied background-error standard deviations are those which would be produced by a Kalman filter implemented over the assimilation window. Although this is one major advantage of 4d-var if the dynamics and simplified physics used in the tangent-linear and adjoint models are accurate, it could be a drawback if this simplified model is not a good approximation of the evolution of the atmosphere. Rabier et al. recommended then to work more and more on the physical parameterization part in the TL/AD, to avoid mismatches between the full nonlinear model and its tangent-linear and adjoint versions.

This recommendation is certainly very important, but efforts must also be done on the regularization of the adiabatic model because in some cases such as the one studied in this paper under some dynamical characteristics and over areas containing no humidity data, model nonlinearities lead to unrealistic increments.

In this study, investigations have been performed on a case of unrealistic specific-humidity analysis over Spain, which happened on April 24th, 2001. The experiments showed that the anomaly is due mainly to AIREP and SYNOP data. AIREPs act through their wind observations and SYNOPs act through their surface-geopotential (pressure) informations. The examination of the terms in the TL/AD formulations responsible for the contribution to these unrealistic increments revealed that AIREP and SYNOP data don't act in the same way, but through two different terms that could be zeroed by imposing dry adiabatic nonlinear and TL/AD models. These terms constitute a great strength of 4d-var allowing multi-variate effects and the use of dynamics as a new source of information to be combined to other ingredients in a data assimilation process. But, they could be very harmful if TL/AD assumption is not satisfie,d especially in presence of high nonlinearities which could be caused by a baroclinic atmosphere evolving rapidly during the assimilation window.

Single-observation experiments enforce this stipulation and indicate in a more comprehensible way the context favouring the occurrence of such phenomena.

The solution of such problem could be, besides the work on enriching data assimilation systems with informations on humidity, to regularize some terms in the nonlinear model in order to avoid the rapid development of nonlinearities, or to add some extra terms parameterizing the second order terms in the TL/AD formulations. This will be the next action to be achieved.

References.

Rabier, F., J.-N. Thépaut and P. Courtier, 1998a : Extended assimilation and forecast experiments with a four-dimensional variational assimilation system. Q. J. R. Meteorol. Soc., 124 , 1861-1887.

Rabier, F., H. Jarvinen, E. Klinker, J.-F. Mahfouf and A. Simmons, 2000 : The ECMWF operational implementation of four-dimensional variational assimilation. I : Experimental results with simplified physics . Q. J. R. Meterol. Soc., 126, 1143-1170.

Rabier, F., and P. Courtier, 1992 : Four-dimensional assimilation in the presence of baroclinic instability. Q. J. R. Meteorol. Soc. , 118, 649-672.

Rabier, F., E. Andersson, J. Haseler, P. Unden, P. Courtier, G. Kelly, D. Vasiljevic, C. Brankovic, C. Cardinali, C. Gaffard, A. Hollingsworth, C. Jakob, P. Janssen, E. Klinker, A. Lanzinger, M. Miller, A. Simmons, B. Strauss, J.-N. Thépaut, P. Viterbo, 1998: The ECMWF implementation of three-dimensional variational assimilation (3D-Var). III: Experimental results. Quart. J. Roy. Meteor. Soc., 124, 1831-1860.

Joly. A., 1992 : ARPEGE/ALADIN : adiabatic model equations and algorithm. ARPEGE/ALADIN Formulation version 0.0.

Rabier, F. and J.-F. Mahfouf, 2000 : The ECMWF operational implementation of four-dimensional variational assimilation. II : Experimental results with improved physics. Q. J. R. Meterol. Soc., 126, 1171-1190.

Klinker, E., F. Rabier, G. Kelly and J.-F. Mahfouf, 2000 : The ECMWF operational implementation of four-dimensional variational assimilation. III : Experimental results and diagnostics with operational configuration. Q. J. R. Meterol. Soc., 126, 1191-1215.