Assimilating remotely sensed cloud optical thickness into a mesoscale model

The Advanced Regional Prediction System, a mesoscale atmospheric model, is applied to simulate the month of June 2006 with a focus on the near surface air temperatures around Paris. To improve the simulated temperatures which show errors up to 10 K during a day on which a cold front passed Paris, a data assimilation procedure to calculate 3-D analysis fields of specific cloud liquid and ice water content is presented. The method is based on the assimilation of observed cloud optical thickness fields into the Advanced Regional Prediction System model and operates on 1-D vertical columns, assuming that the horizontal background error covariance is infinite, i.e. an independent pixel approximation. The rationale behind it is to find vertical profiles of cloud liquid and ice water content that yield the observed cloud optical thickness values and are consistent with the simulated profile. Afterwards, a latent heat adjustment is applied to the temperature in the vertical column. Data from several meteorological stations in the study area are used to verify the model simulations. The results show that the presented assimilation procedure is able to improve the simulated 2 m air temperatures and incoming shortwave radiation significantly during cloudy days. The scheme is able to alter the position of the cloud fields significantly and brings the simulated cloud pattern closer to the observations. As the scheme is rather simple and computationally inexpensive, it is a promising new technique to improve the surface fields of retrospective model simulations for variables that are affected by the position of the clouds. Correspondence to: D. Lauwaet (dirk.lauwaet@vito.be)


Introduction
Mesoscale atmospheric models are used extensively to reconstruct high-resolution regional atmospheric conditions as an input for e.g.hydrological, land surface or air pollution models.Although sophisticated techniques are used to parameterize clouds and precipitation, a large source of uncertainty in the model results remains in predicting the location of cloud systems at high spatial resolutions.As clouds have a strong impact on the surface energy budget and hence the local temperatures, an inaccurate simulation of the overlying cloud cover is problematic for certain applications that need correct surface input data.The assimilation of satellite data into the atmospheric model can play an important role in providing improved model results on a local scale.
Cloud assimilation studies have focused mainly on cloud retrievals from radar data, either with one-dimensional variational schemes (1DVAR) (Benedetti et al., 2003) or with more complex models in 3DVAR (Hu et al., 2006a, b) and 4DVAR (Sun and Crook, 1998;Vukićević et al., 2004).Recently, Benedetti and Janisková (2008) used a 4DVAR system to assimilate Moderate Resolution Imaging Spectroradiometer (MODIS) cloud optical depth observations into the European Centre for Medium range Weather Forecast (ECMWF) model.Their results show a positive impact on certain variables like the distribution of cloud ice water content but the assimilation did not always improve the analysis fit to the observations.However, the large computational infrastructure needed to run and maintain these systems are limiting their use for smaller research centres and universities.
Other simpler and faster methods exist that attempt to retrieve model cloud microphysics from satellite observations or other sources.Soutu et al. (2003)  Their procedure followed the Local Analysis and Prediction System (LAPS, Albers et al., 1996) and clearly improved the model's skill to predict precipitation amounts.Another method is used by Yucel et al. (2003), who applied a nudging assimilation technique to ingest remotely sensed cloud cover and cloud top height data into their mesoscale atmospheric model.The cloud ingestion was found to improve the ability of the model to capture the variation in surface fields associated with cloud cover.However, they suggested that it would be necessary to modify the model dynamics and thermodynamics to be consistent with the ingested cloud fields.
In this context, the goal of the research reported here is to present a rather simple and computational fast cloud assimilation scheme.The scheme applies optimal interpolation with latent heat adjustment for the assimilation of cloud optical thickness (COT) observations into a mesoscale atmospheric model to study the effect on the simulated surface fields associated with cloud cover.The model used for this study is the Advanced Regional Prediction System (ARPS), a non-hydrostatic mesoscale atmospheric model developed at the University of Oklahoma (Xue et al., 2000(Xue et al., , 2001)).Although the ARPS model has its own cloud analysis package (ADAS, Brewster, 1996), it is not used in our study.ADAS needs information on the vertical extent of the clouds to estimate cloud types and related in-cloud vertical velocities which can not be derived from our two-dimensional cloud optical thickness data.
The remainder of the paper is organized as follows.A description of the atmospheric model and a set-up for the model simulations are given in Sect. 2. In Sect.3, the details of the cloud assimilation scheme are presented.The results are evaluated and discussed in Sect. 4 while conclusion are given in Sect. 5.

Numerical model and data description
The Advanced Regional Prediction System (ARPS) includes conservation equations for momentum, heat, mass, water (vapour, liquid and ice), sub-grid scale turbulent kinetic energy and the state-equation of moist air.The finite-difference equations of the model are discretized on an Arakawa Cgrid, employing a terrain following coordinate in the vertical direction.Advection is solved with a fourth-order central differencing scheme and leapfrog time stepping.Turbulence is represented by the 1.5-order turbulent kinetic energy (TKE) model, and the Sun and Chang (1986) parameterization for the convective boundary layer.The 6category water/ice scheme of Lin et al. (1983) accounts for the model microphysics while the Kain-Fritsch cumulus parameterization scheme solves the cumulus convection (Kain and Fritsch, 1990).In order to suppress numerical noise, a fourth-order monotonic computational mixing was applied, following Xue (2000).
Land surface processes are parameterized by the Soil-Vegetation-Atmosphere Transfer model of De Ridder and Schayes (1997).The two primary parameters of the land surface model are the type of vegetation, which is derived from the Coordination Information Environment (CORINE) land cover data and the soil texture, which is assumed to be that of a loamy soil, homogeneous across the domain.Among the secondary parameters, vegetation fraction is based on the normalized difference vegetation index (NDVI) from the SPOT-VEGETATION satellite imagery.
Data with a 0.25 • horizontal resolution from the global operational analysis by the ECMWF are used as initial conditions and as 6-hourly lateral boundary conditions for the model runs.The ARPS model domain has a grid spacing of 16 km and a domain size of 1600 km × 1600 km, centred over Paris (Fig. 1).In all simulations, 35 vertical levels are employed with a grid spacing of 25 m near the surface increasing to 1 km near the upper model boundary, located at 20 km altitude.The simulations are initialized on 1 June 2006 at 00:00 LT and run until 30 June 2006 at 24:00 LT.This month is chosen to test our assimilation scheme as during some periods of this month, the model has problems in simulating the right amount and position of clouds, which leads to large errors in some surface variables as will be shown in Sect. 4.
The cloud optical thickness images for June 2006 are retrieved from visible and shortwave infrared imagery from the Spinning Enhanced Visible and InfraRed Imager (SE-VIRI) onboard the Meteosat Second Generation satellite with a semi-analytical cloud retrieval algorithm (Pandey et al., 2011).This algorithm is based on the retrieval algorithm of Kokhanovsky et al. (2003) for the estimation of cloud optical thickness.The details of the scheme can be found in Pandey et al. (2011).As Meteosat is a geostationary satellite, the algorithm provides COT images every quarter of an hour during daytime (06:00-20:00 LT).These images are assimilated every 15 min into the ARPS model following the procedure that is explained in Sect.3.
To test the effect of our cloud assimilation procedure, 2 m air temperature and specific humidity data for 2 observational stations close to Paris (Melun and Trappes) and a station in Bordeaux are gathered from the National Climatic Data Center (NCDC) dataset (Fig. 1).Furthermore, 2 m air temperature, specific humidity and incoming shortwave radiation data for 3 more stations (Fontainbleau, Grignon and Oensingen, Fig. 1) are taken from the CarboEurope Integrated Project.

Cloud optical thickness assimilation procedure
The data assimilation procedure applied in this study calculates 3-D analysis fields of specific cloud liquid and ice water content (q c and q i ) and operates on 1-D vertical columns.The rationale behind the method is to find vertical profiles  of q c and q i that yield the observed cloud optical thickness fields τ 0 , and that are consistent with the background (simulated) profile, in the sense that clouds are put in layers with a large humidity.This a priori assumption is required as τ 0 does not contain height information.

Background COT
Consider a vertical model column containing n grid cells irregularly spaced at positions z i (i = 1, ..., n), starting from the surface.Each layer (thickness z i ) is characterized by a simulated specific cloud water content q cbi , which can be either liquid or solid (ice) water.Noting that the quantity ρ i q cbi z i is the incremental liquid/ice water path (in kg m −2 ) of layer i (ρ i being the air density of layer i), the incremental optical depth contributed by layer i is given by (Salby, 1996): with ρ w = 1000 kg m −3 the density of liquid water, and r ei the effective radius of cloud droplets in layer i, which is parameterized in ARPS as a function of temperature, to account for the different values of this quantity for liquid versus solid water.The subscript "b" denotes the background fields.A list of all the symbols is provided in Appendix A.
The full model-based columnar optical depth is then given by: H is a so-called observation operator, which linearly maps q cb onto a background (i.e.simulated) optical thickness (τ b = Hq cb ).

Optimal interpolation
Given observations of cloud optical thickness for a given position on the globe, the best linear unbiased estimate of cloud water content is given by (Kalnay, 2003): with q ca the vector containing the analyzed cloud water content values at level i, and q cb likewise containing the values generated by the model ("background" or first guess value).
The gain matrix is given by with B the background error covariance matrix and R (≡σ 2 τ ) the observation error covariance, which in this case is a scalar since τ o itself is a scalar quantity.
We will assume that B is a diagonal matrix.This is not entirely realistic, since errors of cloud water content at different vertical layers may be correlated, especially if these layers are close to each other in comparison to the typical thickness of a cloud layer.Nevertheless, it is difficult to estimate these inter-layer correlations and, moreover, the thicknesses of the layers that are prone to cloud development (sufficiently far away from the surface) are rather thick, thus making this diagonality assumption less of a problem.A diagonal background error covariance matrix has the advantage of leading to a fairly simple final expression for the analyzed specific cloud water content, as shown below.In Appendix B, the results of a test with a non-diagonal B matrix are presented in order to have an estimate of the impact of this assumption.
Thus, in case B is a diagonal matrix with elements σ 2 ci δ ij (with δ ij the Kronecker delta), one has: The main challenge is to specify the σ ci in a suitable manner, in particular in such a way that model layers with a high simulated humidity are more affected than the drier layers.

Cloud water background error variance
The specification of the cloud water background error variance σ ci of each model layer is not straightforward, in particular when a layer contains no simulated liquid or ice water (q cbi = 0).One might be tempted to set σ ci = 0 in such a situation, but from the analysis equation above it is clear that q cai will also be zero then, even if a cloud is observed (τ o > 0).Clearly, a non-zero cloud water background error variance must be assigned, even for non-saturated layers.Simply taking σ ci as a fraction of q cbi will not work for the reasons just explained.The background error variance matrix will therefore be established starting from a probability density function for total specific water content q t , which is defined as the sum of vapour and cloud (liquid/ice) contributions, i.e. q t = q v + q c .It should be noted that by assigning a cloud variance to a non-cloudy background layer, the analysis implicitly allows to adjust the water vapour profile in absence of background clouds.The goal is now to find the cloud water content error, given the error on total water content.The error on the latter needs to be specified a priori, in our case this will be done as a fixed fraction of total water content (see Sect. 3.4).
We employ a normal distribution, given by: with q tb the background (simulated) value of q t , and σ t the standard deviation of the distribution, which is a measure for the error on simulated total water q t .Figure 2 presents the concept of the normal distribution of q t .Cloud water is that portion of q t which is in excess of the saturated value, denoted q s (≡ q sat (T )), with T the temperature of the considered layer, so that q c = (q t − q s )H (q t − q s ), with H (.) the Heaviside step function, which is unity for a positive argument and zero otherwise.Using this, taking the simulated cloud water content q cb as expected value for this quantity, and omitting the layer index i for the moment, the error variance of simulated cloud water can be calculated as follows: such that: with x s = (q s −q tb ) √ 2σ t and erfc(.) the complementary error function.It should be noted that in these formulas the error on the background saturated specific humidity, hence the background temperature field, is ignored to simplify the assimilation scheme.In Appendix C, the effect of this approximation is assessed and it is shown that it does not have a significant effect on the outcome of the assimilation procedure.

Implementation in the ARPS model
In the above, the standard deviation on the simulated total water content and the observed cloud optical thickness are still unknown.These standard deviations are expressed as a fraction of q t and τ 0 , respectively, i.e. σ t =aq t and σ τ =bτ 0. In our study, a value of 0.3 is adopted for coefficient a (i.e.±30 % error of q t ).This value is obtained from a comparison between modelled and observed specific humidity profiles at Trappes, as more than 80 % of the observed data points are within this error margin of the simulated profiles.For coefficient b, a rather conservative value of 0.25 is adopted (i.e.validation study of our cloud optical thickness product with Cloudsat COT data (Pandey et al., 2011).To test the sensitivity of the assimilation procedure to these two coefficients, a sensitivity study is performed by varying the values between 0.1 and 0.5.The results are presented at the end of Sect.4.2.
The resulting q ca from the assimilation procedure is defined as cloud liquid water when the temperature is warmer than −10 • C, and as cloud ice when the temperature is colder than −30 • C. A linear ramp is applied for the temperature in between.Whenever a non-saturated (and cloudless) layer becomes cloudy after the analysis, the specific humidity is set to its saturated value.Last, a latent heat adjustment to temperature based on the added or subtracted amounts of q c and q i is applied, according to the formula: where L v and L f are the latent heat of vaporization and fusion at 0 • C respectively, and C p is the specific heat of dry air at constant pressure.

Results of the assimilation procedure
This section describes the results for a COT assimilation experiment (EXP) for the month of June 2006, compared to a reference simulation (REF) with a setup identical to the cloud assimilation experiment to provide a benchmark for the effect of the introduction of cloud optical thickness data.

Comparison to observations
Figure 3 shows the impact of the COT assimilation on 2 m air temperatures, measured at Melun and Oensingen (locations in Fig. 1).It is apparent that the reference simulation does not correctly capture a sharp temperature decrease halfway through the month (Julian day 166 in Melun and 168 in Oensingen) and keeps on overestimating the temperatures around noon by a few degrees during the rest of the month.This is substantially improved by the COT assimilation which picks up these temperature decreases on the respective days and the following period.The problems for the reference simulation are caused by a wrong location of a cold front and associated overlying cloud cover during these days, as will be shown later on.The assimilation procedure is thus capable of improving the cloud fields and yielding more accurate temperature values.Note that sometimes also the temperature during night time improves in the EXP simulation (e.g. on Julian day 177 in Melun) although the assimilation scheme is only active during day time (as satellite data from a visible channel are needed to derive COT).This is due to the transportation of the assimilated moisture throughout the model domain.
These findings are further demonstrated in Table 1, which shows the results for 4 of the observation stations throughout the model domain.The COT assimilation reduces the rootmean-square error (RMSE) between modelled and observed values for all stations.Also the correlation coefficients between the modelled and observed time series are higher for the assimilation experiment.The positive bias for the stations around Paris that is present in the reference simulation is decreased.However, there is a slight negative effect of the assimilation on the biases in Oensingen en Bordeaux.Overall, the statistics of these time series are clearly improved by the assimilation.
Another variable that is closely linked to the cloud fields, is the surface shortwave radiation.The results for the observation stations of Grignon, Fontainbleau and Oensingen are presented in Fig. 4 and Table 2.As for the temperature, the COT assimilation experiment clearly improves the time series, especially on the problematic Julian days 166 (Grignon) and 173 (Grignon and Oensingen).The statistics show a substantial reduction of the bias and RMSE and higher correlation coefficients.These numbers confirm that the improvement in the assimilation experiment is linked to a better position of the cloud cover in the model.
However, the impact of the COT assimilation is not positive for all variables, as shown in Table 3.The specific humidity at the surface is in good agreement with the observations for the reference simulation, whereas it is overestimated by about 1 g kg −1 for most observation stations when the assimilation scheme is applied.Only in Oensingen, the assimilation improves the underestimated humidity in the reference simulation.The extra moisture is caused by the fact that the assimilation scheme sets the humidity of a layer to its saturated value whenever a non-saturated layer becomes cloudy.As the reference simulation underestimates the amount of clouds compared to the observations, an increase of the hu- midity is the logical consequence in this case.This can not be avoided if we want to retain the assimilated clouds, otherwise they would evaporate instantly.In their cloud assimilation experiments, Benedetti and Janiskova (2008) also  noticed that the humidity was affected in a slightly negative way by the assimilation.In order to evaluate the effect of the assimilation on the vertical distributed moisture, the total water columns from radio-soundings launched at Trappes are compared to the reference and assimilation experiments (Fig. 5).Overall, there is a good agreement between model results and observations, as the mean RMSEs of both the REF and EXP simulations are only 4 mm.The assimilation scheme has only a small impact on the total water column and slightly improves the small negative model bias.So the problems with the humidity that are mentioned before are limited and do not translate in drastic changes in the total vertical moisture amounts.

Impact on temperature, humidity and cloud parameters
To assess the impact of the assimilation procedure on the model variables in the entire model domain, mean zonal differences between the experiment and the reference for temperature, specific humidity and vertical wind speed are shown in Fig. 6.  the temperature and specific humidity.Here, the temperature values of the experiment have a tendency to be lower (up to 0.5 K) while the specific humidity is augmented (up to 0.5 g kg −1 ).Both temperature and moisture changes are in phase to enhance cloud formation.In the upper levels, a slight temperature increase is visible for the assimilation experiment, which is caused by latent heat release during the formation of additional clouds.As shown in the previous section, the temperature decrease near the surface improves the positive bias in the reference simulation, while the moisture increase has a negative effect when compared to the observations.Regarding the vertical wind speed, the largest changes clearly occur in the southern part of the model domain where most extra water is added to the model (see Fig. 7).The enhanced convection and latent heat release causes more updrafts in the higher model layers and downdrafts close to the surface.It should be noted that these changes are rather small (less than 10 %) in comparison to the mean vertical wind speeds.As a response to the changes induced by the assimilation scheme, there is a noticeable redistribution in liquid water path and ice water path in the model domain (Fig. 7).The changes appear to have a rather varied structure over the largest part of the model domain.Most increases occur along the southern boundary of the domain and over the alpine region.Overall, there is a clear tendency of increased cloud amounts in the assimilation experiment.The monthly mean values of the liquid and ice water paths are raised by 25 and 8 % respectively.Regarding the overestimation of incoming shortwave radiation by the reference simulation (Table 2), these changes appear to have a positive impact on the model results.The assimilation also has an impact on the modelled rainfall amounts with a logical positive trend in the regions where more cloud ice is produced (last panel of Fig. 7).The overall effect is relatively large as the monthly mean rainfall is increased by 26 %, although it should be noted that this is a dry month where eventual changes have a strong effect on the overall statistics.
The direct impact of the COT assimilation on the modelled cloud fields is presented in Fig. 8.In this figure the position of the clouds is shown on 15 July at noon, when a cold front passes Paris which is not picked up in the reference simulation (Fig. 3).The top two panels of the figure show a comparison between the cloud optical thickness product of the Moderate Resolution Imaging Spectroradiometer (MODIS) and our product from MSG Seviri.The positions of the clouds in both images are clearly very similar, although the mean COT value of MODIS is 18 % higher than the MSG Seviri value.Finally, the sensitivity of the assimilation scheme to the coefficients a and b, related to the standard deviation of the simulated total water content and the observed cloud optical thickness, respectively, is tested.The results of this sensitivity experiment on the statistics for Melun are presented in Table 4.The values of these parameters have an effect on the model results and the scheme seems a little more sensitive to σ t .Considering the results of both the temperature and humidity at the surface for this station, the initial choice of the values for a and b seems accurate as none of the sensitivity experiments can improve on these results.

Discussion and conclusions
In this paper, a technique has been presented to assimilate cloud optical thickness information into a mesoscale atmospheric model to yield improved diagnoses of surface solar radiation and temperature.The technique comprises an optimal estimation of cloud liquid and ice water in 1-D vertical columns together with a latent heat adjustment.The scheme requires some assumptions including an independent pixel approximation, but it is rather simple and computationally inexpensive, especially when compared to the 4DVAR systems that are currently developed (e.g.Benedetti and Janisková, 2008).The goal of the assimilation scheme is to improve retrospective model simulations by feeding the model with observed cloud optical thickness images every 15 min.
Results for the month of June 2006 show a positive impact of the assimilation on near-surface temperatures and incoming shortwave radiation, two variables that are closely linked to the overlying cloud cover and are crucial as input in, for instance, air pollution models.However, comparison to specific humidity observations show that the changes induced by the assimilation do not always improve the model fit to the observations.The assimilation scheme tends to induce overestimations of humidity close to the surface due to the fact that a layer is set to saturation when it becomes cloudy.This is necessary to retain the new clouds in the model and the same technique is used in the cloud analysis scheme of Soutu et al. (2003).Although the moisture field in the lowest 2000 m of the model domain is affected in a slightly negative way, the results show that the position of the cloud fields are more accurately simulated when the cloud observations are assimilated.
We can thus conclude that it is feasible to introduce the presented COT assimilation procedure in a mesoscale atmospheric model.The results show that the benefits of the assimilation to the surface temperature and radiation fields outweigh eventual inconsistencies that are caused by only adjusting the moisture and temperature fields of the model.As the procedure is simple and fast, it is a promising new technique to improve the quality of surface level model output of retrospective simulations.transposed observation operator matrix K gain matrix L f latent heat of fusion (J kg −1 ) L v latent heat of vaporization (J kg −1 ) q specific humidity (kg kg −1 ) q c specific cloud liquid water content (kg kg −1 ) q ca analyzed cloud water content (kg kg −1 ) q cb simulated (background) cloud water content (kg kg −1 ) q i specific cloud ice water content (kg kg −1 ) q s saturated specific humidity (kg kg −1 ) q sb simulated (background) saturated specific humidity (kg kg −1 ) validation data.The ARPS meteorological model was made available by the Center for Analysis and Prediction of Storms at the University of Oklahoma.We would also like to acknowledge the European Environment Agency for making available the CORINE land-cover data, the US Geological Survey for the GTOPO30 dataset, the Jet Propulsion Laboratory for the SST imagery, NASA's Earth Science Division for the GLDAS data, and our VITO colleagues at TAP for the SPOT-VEGETATION NDVI imagery.Meteorological forcing data were obtained from the European Centre for Medium-Range Weather Forecasting (ECMWF).

List of symbols
Edited by: A. Geer constructed cloud fields for their forecasts over the Galician Region in Spain based Published by Copernicus Publications on behalf of the European Geosciences Union.D. Lauwaet et al.: Assimilating cloud optical thickness on relative humidity values from the NCEP Aviation Model.

Fig. 1 .
Fig. 1.Location of the model domain and the observational stations.

Figure 2 .Fig. 2 .
Figure 2. Probability distribution function for total water content (blue line), with expectation 2 value qtb, and standard deviation σt.The light shading corresponds to the area above the 3 saturated specific humidity (denoted qs), which contains cloud water.4 5 6 7 σ

Fig. 5 .
Fig. 5.The water column at Trappes from radio-sounding measurements (dark grey bars), the Reference simulation (black bars) and the COT assimilation experiment (light grey bars).

Figure 6 :
Figure 6: Monthly mean differences between the COT assimilation experiment and the Reference simulation, averaged per latitude band.Fig. 6.Monthly mean differences between the COT assimilation experiment and the Reference simulation, averaged per latitude band.

Fig. 7 .Fig. 8 .
Fig. 7. Monthly mean differences in liquid water path (upper panel), ice water path (middle panel) and rainfall amounts (lower panel) between the COT assimilation experiment and the Reference simulation.

a
coefficient in relation to σ t (−) b coefficient in relation to σ τ (−) B background error covariance matrix C p specific heat of air at constant pressure (J kg −1 K −1 ) erfc complementary error function H observation operator matrix H heaviside step function H T

Table 1 .
Statistics of the 2 m air temperature for the entire month of June 2006.

Table 2 .
Statistics of the incoming shortwave radiation for the entire month of June 2006.

Table 3 .
Statistics of the 2 m specific humidity for the entire month of June 2006.

Table 4 .
Statistics for Melun for all the sensitivity experiments.For the reference simulation, the clouds are positioned too far to the east and there is no strong front structure visible.The COT assimilation scheme is able to alter the cloud structure significantly and brings it much closer to the observations over the central region of the model domain.Although the scheme is not able to get rid entirely of the overestimation of clouds in the west side of the model domain, it is clearly able to improve the model simulation.