"Mesoscale variational assimilation for land surface variables"

Gianpaolo BALSAMO : rapport for the Newsletters 20&3, November 2001


A variational assimilation technique is considered to analyse the soil variables by assimilating 2m data of temperature and relative humidity. We consider 2-D (z and t) variational approach, under the assumption of horizontal decoupling of surface processes between grid points. This hypothesis is firstly supposed considering the small scale involved in the soil processes far from the actual grid mesh in NWP (few meters against few km). A validation test is then designed to evaluated the influence of neighbours grid points. Thus, the method is applied on every grid point singularly and the gain matrix is directly computed given the small dimension of the problem. The variational technique keeps count of the full physics of the model and therefore the corrections applied on the control variable are adapted to the current meteorological conditions and the grid point characteristics (texture, vegetation fraction, LAI, ... and the previous soil state). The observation departures are computed at the synoptic time (00 06 12 18 UTC) and a linear estimation of observation operator is obtained by applying a perturbation to the control variable (Ws, Ts, Wp, Tp). A first test on a 48 h time window for the total soil moisture content (Wp) analysis is presented. A sequential assimilation cycle every 6h is also studied as it allows to focus on specific problems.

The linearized variational approach

The assimilation of the data in NWP consists in best fit the data and the model previous information. This optimal state is obtained by the minimization of the cost function J(x), of the form:

J(x) = J b(x) + J o(x) = ½ (x - x b)T B-1 (x - x b) + ½ (y - H(x))T R-1 (y - H(x)) (1)

where:
x is the control variables vector (x = (xi) 1 £ i £ n )
x b is the control variables vector (x b = (xib) 1 £ i £ n )
y is the observation vector (y = (yi) 1 £ i £ n )
H is the observation operator, H its linear approximation
x t is the real state, y t = H(x t )
B is the matrix of forecast errors covariances
R is the matrix of observation errors covariances

GB1.gif

The observation operator H transports the state vector x into the observations space. This includes geometrical and physical interpolation as well as integration in time. This method needs normally an effort of coding the tangent linear and the adjoint of the direct code of the model in order to perform the minimization of the cost function. For low dimension problems, the minimization of the cost function can be imposed directly by the expression:
Ñx J(x) = 0 = B-1 (x - x b) + HT R-1 (y - H(x)) (2)

When this condition is satisfied x = x a . The H operator can be expressed by the Taylor prime order expansion under the linear hypothesis:
H(x+dx) = H(x) + Hdx (3)

By substituting the expression (3) in the equation (2) and rearranging it, is obtained:
Ñx J(x a) = 0 = B-1 (x a - x b) + HT R-1 (y - H(x b) - H(x a - x b)) (4)
x a = x b + BHT(HBHT + R)-1(y - H(x b)) (5)

where the terms:
K = BHT(HBHT + R)-1 y - H(x b) (6)
are respectively the gain matrix K and the innovation vector. This equation recalls the formulation of the Best Linear Unbiased Estimate (BLUE) but including the integration in time. The method requires an extra integration of the model initialised with a perturbed state of the control soil variable to analyse. The perturbation provides a linear estimation of the Jacobian of the observation operator H. The linear hypothesis is best satisfied in the vicinity of x b, so the method is suitable to be applied for small correction of the first guess.

Application of the method for Wp analysis

The analysis of total soil water content Wp of the ISBA scheme using 2m observations (temperature T and relative humidity Hu) is presented as an example. Considering the 2m variations between the guess G and its perturbation G':GB.gif

The matrices B, R, H, used for the direct computation of K are simply:

GB2.gif

The innovation vector is given by the difference between the observation and the reference at the grid point for T2m and Hu2m. Here no approximation is required, but the difference is performed at the grid point instead of the observation point. This practically implies a previously analysed field for T2m and Hu2m. The form of the innovation vector is as follows:

GB3.gif

The experimental frame to test the convergence

To investigate the general convergence of the method, several experiments are performed. A model run initialised with a well known reference state of soil moisture Wref is used to extract simulated observations (OBS=GUESS of the reference run). The error of the analysis is evaluated as:

GB4.gif (8)

The accuracy of the method to retrieve an error equal to the perturbation applied is fair. For a perturbation of dWp = 10 mm applied on each grid point, the use of simulated observations allows to retrieve a correction of 10 mm with a mean absolute error lower than 5.2 10- 6 mm. The accuracy is dependent on the error covariances assumed for the background and the observation terms. Therefore to have a greater convergence towards the simulated observations, the standard deviations (s) for observational errors are chosen of low magnitude (0.01 K for T2m and 0.1 % for Hu2m), and that for background (Wp) is put to 100 mm. A set of randomly generated state of Wp with given mean and standard deviation is now considered as first guess. From the different departures (W1, W2, W3 for wet guess, D1, D2, D3 for dry guess) generated with a linear random distribution around a given mean value it is checked the convergence towards the reference (Wref ).

The perturbation G' is produced from the guess G according to the following condition:

GB5.gif

Coupling files are taken from the 16-17 June 2000 and the observations are generated from the reference run without noise (perfect observations).

Assimilation on 48 hours time window

A 48 hours test from the different departures gives indication of mean convergence. The maps of Wp errors show some patterns of not convergent points over eastern Europe and southern Italy that can be put in relation with the synoptic situation. The main reasons for the poor local convergence can be attributed to the little information content in meteorological perturbed areas or in mountainous regions. Anyway, by iterating the convergence it is reached an equilibrium state (from every Wp departure) where the percentage of total land points that converge with an error lower than 10 mm is around 70 % and becomes greater that 94 % if we consider errors lower than 30 mm.

A 6h sequential assimilation cycle

Following the results of the 2-D (z and t) variational approach for a long time window, the method has been applied in sequential way with a 6 hour forecast cycle. This allows to better understand the relations with meteorological situation and to see the impact of each network assimilated. The H matrix is evaluated considering 6h forecast (at time t) starting from the guess G and its perturbation G' (initialised with the corresponding values of Wp at the time t-6h). The convergence at the end of the 48 hours period (after 8 cycles) is reached with an error lower than 10 mm for 82 % of the points and 96 % if we consider errors lower than 30 mm.

Study of the behaviour of the H matrix

The estimation of H matrix by mean of a perturbed forecast could lead to undesired variation of other near-surface or vertically cumulated fields. Particularly the precipitation (as shown in figure 1), the 10m wind, and the incoming solar radiation (and also the cloudiness to cover nocturnal hours) are considered. A further test is done to put in evidence the effect of the synoptic situation on 2m temperature during the integration. The figure 1 shows the time-step evolution of the guess and the perturbed run for T2m at 2 points located respectively in clear sky and rainy situation. Although the mean behaviour of T2m is similar, with a higher temperature for perturbed guess field, the resulting gradient is of opposite sign due to stronger oscillation in rainy conditions. This leads over certain areas (figure 2 left) to a noisy representation of the H matrix, which is responsible for the local poor convergence. Those situations are likely to occur when the atmosphere is interested by precipitation events, reduced solar radiation flux (cloudy condition) or strong lateral advection. A set of constraints are applied by masking grid points where the sensitivity to initial conditions (evaluated as difference between the guess and its perturbation G - G') overpass a given threshold. As a result the main noisy patterns shown in figure 2 (right) are masked out.

Figure 1: 2m temperature evolution in clear sky (Western point :47.27;3.10, left curves) and rainy situation (Eastern point: 53.70;17.01, right curves)

Figure 2: H matrix term at 12UTC of 20000616 - DT2m / DWp (K / m); left: no mask; right: mask of sensitive zones

Validation of the horizontal decoupling hypothesis

A validation test is designed to prove the horizontal decoupling hypothesis. The reference state of soil moisture Wref is taken as the medium value SWI=0.5 over the entire domain of ALADIN-France and simulated observations are generated with a 6 hour forecast. The guess (with a modified soil moisture Wmod) is taken as SWI=0.75 within a box of about 400 ´ 400 km centred over France and SWI=0.5 elsewhere. Inside this main box four other boxes are placed (size 10 ´ 10 km at SWI=0.25, 20 ´ 20 km at SWI=0.25, 40 ´ 40 km at SWI=0.5 -equal to reference-, 60 ´ 60 km at SWI=0.25), to test different scales and to observe the influence on grid points close to the boxes. The model is initialised at 6 UTC and 6 hour forecast is produced. Results are plot in figure 3. The initial difference between the guess and the reference for the soil moisture is plot in figure (a). The difference between the guess and the reference in the 2m temperature and relative humidity at 12 UTC are plotted in figures (b) and (c). The analysis of soil moisture is then obtained at 12 UTC with the constraints elaborated.

Figure 3: Test to assess the horizontal decoupling at 20000616 at 12 UTC for a modified initialisation of Wp at 6 UTC. The initial span Wmod - Wref is shown in (a). (b) and (c) report the DT2m and DHu2m obtained after 6h integration (innovation vector for the soil moisture analysis). The residual error after the analysis done at 12 UTC (Wanalysis - Wref) is given in (d)

The results show a net impact of the structure of Wp on the T2m and Hu2m fields with conservation of all scales features both for T2m and Hu2m. The boundaries present some extended signal of lower magnitude. The intensity of the signal depending on the horizontal scale of the box confirms also some small influence. Finally the analysis is able to retrieve the reference with good approximation in the inner box and reduced overshooting effect on the boundaries.

Conclusions and perspectives

It has been evaluated the use of a linearized variational technique for the total soil water content analysis from 2m simulated observation and for different time windows. The results are encouraging and further tests are necessary to validate the method with real observations. Particularly interesting is the 6h time window as it allows comparison with the statistical OI method. The study of the H matrix highlights several features related to synoptic situation. The mask of sensitive zones allows to further improve the convergence towards the reference state. In future work, the dynamically derived coefficients of correction provided by the variational method will be compared to the statistical OI coefficients used operationally.




Home