Abstract
This study demonstrates successful variational retrieval of land surface states by assimilating screen level atmospheric measurements of specific humidity and air temperature. To this end, the land surface scheme is first validated against the Oklahoma Atmospheric Surface Layer Instrumentation System measurements with necessary refinements to the forward model implemented. The retrieval scheme involves a 1D land surface-atmosphere model, the corresponding adjoint codes, and a cost function that measures residuals between observed and modeled screen level atmospheric temperature and specific humidity. The retrieval scheme is robust when subjected to observational errors with magnitudes comparable to instrument accuracy and for initial guess errors larger than typical model forecast errors. Using varying assimilation window lengths centered on different periods of a day, the sampling strategy is assessed. The daytime observations are more informative compared to nocturnal observations. An assimilation window as narrow as four hours, if centered on local noon, contains comparable information to an expanded window covering the whole day. There exists an optimal assimilation window length resulting from the contest between degrading forecast accuracy and increasing information content. For an assimilation window less than two days, the “optimal” assimilation window length is inversely proportional to the data ingesting frequency.
1. Introduction
Variational data assimilation of energy fluxes is a core component of numerical models used to investigate convective initiation. For example, the hurricane landfall was found (personal communication with L. Leslie, 2014) not only to be sensitive to SSTs but also to be the land surface condition in the coastal region for both flat and hilly terrains. A major obstacle to improving numerical modeling of near surface variables is the sensitivity of model predictions to the accuracy of specification of key land surface parameters. The motivation for retrieving initial land surface conditions and the reason for choosing a variational approach have been justified in [1], where we also discussed the assimilation of ground surface temperature for retrieving initial land surface condition. In reality, except for the data rich areas such as Oklahoma, the only source for skin temperature is that derived from satellite remote sensing. Under cloudy conditions, there is no guarantee for high quality data availability if thermal infrared sensors are used. If microwave sensors are involved as the sources for skin temperature observations, the accuracy of the measurements remains an ongoing problem and there only will be a couple of satellite overpasses per day from one low Earth orbit platform [2]. Data collected this way may be too sparse for many study purposes. It would be desirable if routinely available screen level observations can be utilized to infer initial land surface conditions.
Here we discuss several closely related numerical modeling studies of recent years so as to put our work into the appropriate context. Indirect methods have been proposed to retrieve information on soil moisture from atmospheric information and conventional land use data. Mahfouf [3] attempted the retrieval of surface moisture from observations of screen level air temperature and relative humidity (10 m wind and 2 m temperature and humidity as SYNOPs). Using the Interactions Soil Biosphere Atmosphere (ISBA [4]; NP89 henceforth) land surface scheme, positive conclusions were obtained for an agricultural crop area under clear sky conditions. Two possible approaches were described: a variational algorithm where a cost function is minimized over an assimilation period and a sequential assimilation scheme that consists of a set of predictions and in-static correlations of soil moisture. Both methods were validated against in situ data collected during the HAPEX-MOBILHY [5] field experiment, using a one-column model to represent the interactions between surface processes and the planetary boundary layer structure. Bouttier et al. [6, 7] further argued for the feasibility of estimating both superficial and deep/bulk layer soil moisture using the time evolution of atmospheric temperature and relative humidity near the surface. They also emphasized the fact that this method requires a close relationship between the near surface atmospheric state and the soil moisture.
Calvet et al. [10] also applied a nonvariational inverse scheme for estimation of the bulk soil moisture content using surface variables (either surface soil moisture or temperature). They argued that knowing the atmospheric forcing (especially precipitation) and four to five surface soil moisture measurements over a two-week period were adequate for retrieving the bulk soil moisture by inverting ISBA scheme. Because there is a strong relationship between the deeper layer soil moisture and surface soil moisture, especially when the vegetation cover is in full growth, it is feasible to infer the bulk soil moisture by minimizing the error in the predicted surface soil moisture. For assimilating surface temperatures, they realized that satisfactory retrieval is obtainable for only relatively dry conditions.
Based on the same forward model as used by [10], Bouyssel et al. [8] performed a series of variational surface analyses using screen level atmospheric parameters. Their tangent linear analysis identified several interesting topics that inspired the present study. Their study, however, was based solely on synthetic data (i.e., identical twin experiments).
Compared with Ren [11], the new set of experiments described here is more indirect and more demanding because screen level atmospheric measurements are assimilated to infer the initial soil model conditions [3, 8, 12]. In addition to synthetic data, the real Oklahoma Atmospheric Surface Layer Instrumentation System (OASIS) observations also will be assimilated.
Since the data set used and the land surface component of the forward model are identical to those of [11], our description focuses on describing the atmospheric planetary boundary layer scheme and the forward model verification followed by a rather general description of the variational formalism and the construction of the adjoint model. The retrieval experiments assuming data availability of screen level atmospheric parameters are performed and results are analyzed for both synthetic observations and real OASIS measurements.
2. Data and the Forward Model
There are two prerequisites for a successful variational retrieval: an accurate forward model (in terms of physical mechanisms, numerical schemes, and coding) and high quality observational data and a cost function sensitive to the control variables. If the uncertainty in a retrieved variable contributes little to the model-data misfit, as measured by the cost function, that variable simply cannot be effectively retrieved. Thus, effort must first be applied to verify the forward model and rectify any problems related to soil temperature/moisture estimation. During the forward model calibration and backward model inversion, the high quality, high-frequency OASIS measurements are used as the ground truth and as input for retrieval experiments. The second requirement as applied to the retrieval of soil moisture and temperature from atmospheric forecasting errors is that the impact of soil moisture/temperature on near surface observations dominates the impacts of other error sources [6, 7].
2.1. Forward Model Description
The forward model is based on the Advanced Regional Prediction System (ARPS [13–15]). All of our experiments, including forward prediction and 4D-Var retrieval, will be conducted in 1D column mode. This is necessary because, otherwise, a full 3D adjoint code would be required for the atmospheric component of model, and a full blown 4D-Var data assimilation system is involved. This is a major undertaking that is beyond the scope of this study primarily because the associated computational cost would preclude routine implementation.
Almost all land surface models used in atmospheric prediction systems are column based; namely, they do not include horizontal transport of heat or moisture. This is also true of the land surface model in the ARPS that we use. For our study, the vertical boundary layer mixing process is of foremost importance in the atmospheric model, particularly, for the periods when weather is inactive so that horizontal advection is less important. In fact, for the experiments that use observational data (real data experiments), the quiescent atmospheric periods are intentionally chosen, so that our 1D assumption is valid.
The core of the forward model involves solving the coupled energy and water budget equations of the land surface and of the overlying atmospheric boundary layer (Figure 1). Because of the 1D assumption, the hosting system (ARPS) is trimmed to keep a minimum number of relevant components. In all prognostic equations of the atmospheric model, the advection terms are neglected and so are horizontal mixing terms. For momentum equations, the Coriolis terms are retained, and the horizontal pressure gradient terms are expressed in terms of geostrophic winds. For the potential temperature equation and the water vapor equation, only the vertical mixing term is kept. Within the atmospheric boundary layer, the turbulent eddy coefficients for momentum, heat, and moisture are diagnostically prescribed using the PBL model of Hong and Pan [16]. Since the land surface component has been discussed in [11], our discussion focuses mainly on the atmospheric component.
The PBL is the layer of atmosphere that directly interacts with the land surface. For the retrieval of soil state variables using surface atmospheric observations, the soil component will be coupled with the atmospheric model and the PBL mixing is the most important process. The governing equation readswhere is time, is altitude, is the zonal velocity component (eastward positive), is the meridional velocity component (northward positive), is vertical velocity, and is virtual potential temperature, which is defined as ; here is the potential temperature for dry air, and is the mixing ratio of water vapor. Here is the eastward component of the geostrophic wind, is the northward component of the geostrophic wind, is the Coriolis parameter, is the radiative cooling/heating rate, and is the sink term due to condensation (varnishes for clear boundary layer). Here the overbar “—” denotes the grid resolvable or mean quantity, while the superscription prime “′” indicates subgrid quantity or perturbation. The vertical turbulent flux of a quantity (, , or momentum) is expressed as the covariance between its perturbation and the vertical velocity perturbation. Specifically, and are components of the Reynolds stress tensor, respectively, in the directions east and north; and and are, respectively, the components of the turbulent heat and moisture fluxes.
The Reynolds stress term is usually parameterized as proportional to the vertical gradient of the mean flow to close the system. The key to a turbulent closure is the determination of the mixing coefficient which is usually parameterized using PBL parameterization scheme inside the convective boundary layer and is estimated using a subgrid scale turbulence scheme above the convective boundary layer. In this study, for simplicity, is a diagnostically obtained profile following the work of Hong and Pan [16]. The nonlocal PBL parameterization scheme implemented into ARPS is described in the Appendix.
For our 1D model run, for simplicity, the geostrophic winds for atmospheric levels above the PBL are interpolated from the two subsequent soundings about 12 hours apart and the winds are linearly interpolated to the corresponding model vertical levels. For the layers within the PBL, a constant geostrophic wind profile is assumed, with the wind speed values at the top of the boundary layer. The radiation heating/cooling processes and the microphysics are retained as in ARPS.
Implementation of this PBL scheme in the ARPS framework is not straightforward because the ARPS equations are first written in a Cartesian coordinate projected onto a plane tangent to or intercepting the Earth’s surface. A coordinate transformation into a curvilinear coordinate system is then performed to put the governing equations into an equally spaced computational domain. The special curvilinear coordinate system that ARPS applies isEquation (2) represents a transformation that maps a domain with a vertically stretched grid and an irregular lower boundary to a regular rectangular domain with equal grid spacing in each direction. We call the latter the computational domain. The dynamic equations are discretized in the computational space. Our column run voids the horizontal coordinate transformation and makes the map projection transformation irrelevant or alternatively makes the map projection factor equal to unity. With the definition of , and , the governing equations for atmospheric components becomewhere is the speed of sound, is the ratio of the gas constants for dry air and water vapor, and is the radiative forcing.
In ARPS, wind components and the state variables are defined as the sums of the base-state variables (those with an overbar in (3a)–(3e)) and the deviations from the base state (with primes removed). The base state is typically constructed using an external sounding and is assumed to be horizontally homogeneous, time invariant, and hydrostatically balanced.
is fixed once the elevation and vertical grid setting is set. For all the following numerical experiments, 80 grid points are used in the vertical direction. The vertical grid is stretched from 4 m at the bottom to 396 m at the top level, according to the hyperbolic tangent function given by Eq. in [13]. The vertical dimension of the simulated domain is ~32 km, deep enough for assuming zero pressure perturbation at the top of atmosphere. This is actually the upper level boundary condition for numerically solving the hydrostatic perturbation pressure equation (3e).
For a description of the fully fledged 3D ARPS grid stencil setting, the reader is referred to [14, 15]. For our discussion, we are concerned only about its vertical grid setting (see left side of Figure 1). It is a reduced Arakawa C-grid [17], where all prognostic scalar variables are defined at the center of the grid box, while the normal velocity components are defined on their respective box surfaces. Other derived variables are evaluated at locations that minimize spatial averaging in the difference operations. According to the variable arrangement relative to the physical boundary, the second level of scalar variables is the first level above ground surface. The vertical gradient of a quantity is evaluated using the first and second layer and is centered on the land surface, where the vertical turbulent heat flux, H3, is defined.
Since our column does not support gravity waves, the pressure perturbation equation is further reduced to the hydrostatic equilibrium form. The separation of acoustically active and inactive modes as appears in the split-explicit time integration scheme of 3D ARPS is unnecessary. The horizontal processes (i.e., the Coriolis’ force and the pressure gradient force) can be treated explicitly without stability constraints. The leap-frog scheme is used for the time integration. However, for the often large vertical mixing coefficients used by the PBL scheme, vertical turbulent mixing often results in a linear stability constraint on the time step size when treated explicitly, especially when the vertical resolution is high. To overcome this potentially severe restriction on the time step size, the implicit Crank-Nicolson type scheme is used for the vertical mixing so that the time integration is absolutely stable for the mixing terms [14].
In ARPS, the atmospheric and land surface components use different conventions. For example, the equation for potential temperature (see (3c)) is not in the energy form (the heat capacity factor is not shown). Hence, the sensible heat flux output from the land surface model is not of the same quantity as the potential temperature flux at the lower boundary. Similarly, the equation for specific humidity (see (3d)) is written in mass flux form (vapor mass per unit time per unit area), whereas the land surface model provides latent heat flux (energy per unit time per unit area) at the lower boundary. Due to these specific features of the atmospheric model, proper adjustments must be made before surface fluxes provided by land surface model may be used by the atmospheric component of the model. In addition, the formulation of the land surface model uses differences in temperature, while the eddy turbulent mixing term is expressed in potential temperature form; a transform factor should also be taken into account, so that Here is surface layer exchange coefficient, is the horizontal velocity vector for the lowest model layer, is the latent heat of evaporation, and subscripts “sfc” and “” identify surface and the first atmospheric model layer quantities, respectively.
The kinematic surface fluxes are given bywhere and , respectively, are zonal and meridional components of the vector wind . Note also that all fluxes are defined as positive when directed upward in the atmospheric model, whereas they are defined positive when contributing to the ground’s substrate.
In ARPS, the parameterization scheme for estimating and together with the original TKE scheme for estimating in the PBL layer and the force-restore type of soil scheme was tested as in the framework of ARPS against several field experiment data sets. We follow the ARPS estimation of and rather than using those in the original Medium Range Forecast (MRF [18], Appendix) code because we believe the consistency of the original soil scheme with the surface layer parameterization is more important than the consistency with profile in PBL layer.
ARPS diagnoses the PBL depth based on the virtual potential temperature profile. The top of the PBL is assumed to be at the level where the environmental virtual potential temperature exceeds that of the first level above ground. If the atmosphere is stable right above ground, the PBL depth is set to the thickness of the layer below the first scalar point above ground. In this study, the iteratively obtained in MRF is used instead, which includes further adjustment to the height of the inversion by thermal and nonlocal processes.
2.2. Verification of Forward Model
Because the forward model is used as a strong constraint, its accuracy affects the accuracy of the retrieval. This section contributes to verifying the forward coupled land surface atmospheric model. For the land surface quantities (i.e., soil temperature, moisture, LE, and H) the model output is checked against corresponding OASIS measurements. Atmospheric profiles of potential temperature and specific humidity are compared with soundings that are available every 12 hours. Our verification experiments are performed for a variety of situations and the desirable results are obtained. In the following, we present results from one clear sky dry period, August 6–8, 2000. Model forecasts continuously for two days, with a time step of 1 minute for both atmospheric and soil components. The atmospheric model parameters used in this study are listed in Table 1.
The atmospheric component is initialized using sounding of 12 Z on 6 August, 2000. The surface and deep soil moisture values are initialized using the measurements at 5 and 25 cm, respectively. For the selected period, the respective values for superficial and deep soil moisture are 0.253 and 0.278 m3 m−3, respectively. The surface temperature is initialized using a measurement by an infrared instrument (295.8 K), while the deep soil temperature (298.5 K) is initialized using the preprocessed value according the procedure described in [9].
With the above parameter settings, satisfactory descriptions of the soil temperatures are obtained with our model (Figure 2). The surface temperature is predicted accurately, with maximum errors generally less than 2 K. There lacks apparent phase error and a slight warm bias is present in the prediction that is less than 1 K. The time series for deep soil temperature () indicate that revisions to the force-restore temperature equations are very necessary. Otherwise, the deep soil temperature will drift upward and finally assumes the same daily average temperature as that of the surface temperature. Particularly satisfying are the accurate predictions of the surface and deep soil temperatures: the maximum model-data differences of 2 K and 0.5 K, respectively, fall well within the measurement error range. Without the modifications introduced to the model by Ren and Xue [9], the errors would be much larger.
The out-of-phase behavior in the simulated surface soil moisture with respect to diurnal atmospheric forcing is believed to be due to more complicated vegetation activity and the hydraulic displacement of the soil water potential [19–21]. It is beyond the capacity of the current force-restore model, which does not include the effects of hydraulic lift. Without implementing the effects of hydraulic lift, the simulated surface soil moisture shows phase error compared with the observations. However, the absolute differences between forecasts and observations are generally small, typically the magnitude of the measurement error. We emphasize the fact that the amplitude of the daily cycle is rather accurately modeled. The simulation of the deep soil moisture is accurate.
The measured and predicted specific humidity is not exactly the same from the starting time, due to the fact that the model is initialized using a 12 Z sounding and the 2 m value is interpolated from the sounding which is ~1 g/kg larger than the corresponding measurement. The evolution over time occasionally deviates from the measurements. The downward wiggles (15 Z and 21 Z) corresponding to switching on and off of the stability measured by the bulk Richardson number. Using a value of bulk Richardson number < 0.5 helps reduce these wiggles. It is interesting that the measured specific humidity increases during heating period and decreases during nighttime. The out-of-phase features are connected with the surface soil moisture simulation. The absolute difference is generally less than 3 g/kg and the simulation is satisfactory on daily basis. Whether or not the out-of-phase feature affects our retrieval based on OASIS screen level observations depends critically on the sensitivity of specific humidity to soil moisture content. The retrieval scheme can perform at best to reduce the initial guess errors in initial soil moisture values within the perturbation range which can cause specific humidity bifurcating as large as 3 g/kg for this relatively dry case.
Figure 3 indicates that the simulations of LE and are accurate, with peak value differences less than 50 W m−2 during daytime heating period. During nighttime, the model overestimates sensible heat flux, while it underestimates latent heat flux. However, both magnitudes are small (less than 15 W m−2) and well within the instrumental error ranges for such variable measurements.
Accurate prediction of the atmospheric state is also critical for our retrieval experiments. We compare the vertical profiles of atmospheric potential temperature and specific humidity with the corresponding soundings (Figures 4(a) and 4(b)). 12 hours later (00 Z August 7), the modeled and observed potential temperature profiles are close to each other (modeled is slightly warmer) and the difference is generally less than 2 K. 24 hours later, the profiles are still quite similar as to the PBL height, whereas the difference can be as large as 5 K. Most of the errors in the specific humidity profile are near the surface, whereas the prediction is pretty accurate at higher levels. Although the detailed structures in the profile are not predicted after 24 hours, the general trend is quite accurate (<2 g/kg).
(a)
(b)
The time-height cross sections as plotted in Figure 5 illustrate the daily cycle of eddy diffusivities coefficients (for vapor and thermal) calculated using the nonlocal scheme as implemented in our model. It can be seen that the diffusivity coefficients gradually increase from sunrise up to mid-afternoon and then decrease relatively quickly near sunset. The location of the maximum values of and is about 1/3 of the overshoot height. The shapes of ’s profiles experience a transition from like one-side distribution to biased two-sided distribution and to normal distribution as vertical stratifications transit from stable to neutral and to convective. Such results are in concert with those of Hong and Pan [16] (see Figure 2 therein).
A sequence of hourly plots (not shown) indicates that a mixed boundary layer is fully developed by 17 Z. Due to the vegetated ground surface, the superadiabatic layer near the surface is not apparent, whereas the negative potential temperature perturbation develops at the top of the boundary layer as a result of entrainment. The depth of the planetary boundary layer is accurately predicted. The nighttime very cold and dry spikes which appear for the potential temperature and mixing ratio simulations, respectively, are successfully eliminated by implementing the flux distribution scheme of ARPS, which distribute the estimated surface radiative cooling caused near surface air potential temperature drops during nighttime into a certain depth. The depth for swapping flux redistribution was set to 200 m, which is clearly shown in the 3 Z profiles. We are pretty sure that the development and collapse of the planetary boundary layer are accurately described by the model because, on 12 Z August 8 (the ending time of forward integration), the soundings for mixing ratio and potential temperature are rather close to our modeled profiles.
The accurate description of daytime growth and nighttime collapse of the PBL as well as the evolution of the PBL height is important not only because PBL height exhibits a strong daily cycle driven by the surface heating but also because it grows in conjunction with the entrainment of the warm (higher in potential temperature) and dry air from the overlying free atmosphere. The entrainment fluxes, which serve as the main connection between the coupled land surface PBL system and the overlying free atmosphere, warm and dry the mixed layer. Considering that the similarity theory is applied on the mixed layer property and the corresponding surface property, this process of entrainment directly affects not only the PBL energy and moisture distribution but also ultimately the surface energy and moisture budget.
The results confirm the accuracy of the forecast prediction model, which is essential for the retrieval experiments to be conducted based on this forward model and its adjoint. Also, the reader is cautioned that, due to technical difficulties in adjoint coding, some parameterizations in the PBL component as well as in land surface component will be altered. However, all such modifications are acceptable only when they do not alter the overall forecasting performance as shown in this section.
3. Definition of the Cost Function and Verification of the Optimization System
For the soil state variable retrieval problem using screen level atmospheric observations as the constraint, the cost function is constructed as a quadratic function of screen level atmospheric forcing variables. They are related through model dynamics to the land surface model variables. We define , and where is the number of observations during an assimilation period; represents temperature (superscript means matrix transpose), with subscripts “air” indicating the screen level of air. Here is water vapor specific humidity, with the notational conventions as those for temperature. ’s represent the relative confidence accredited to each observation and prediction pair, which are typically the standard deviation of the observational error. is the prior estimate of the initial condition vector and is the prior error covariance. The three ’s weight the model-measurement misfit and prior estimate misfit, respectively, based on the uncertainty involved in each component. For example, if the measurement error was very small and the prior knowledge of the initial conditions was poor, then we would expect the estimation procedure to change greatly the initial guess or prior estimate in an effort to fit the model output closer to the measurements.
When many high quality observations are available in the assimilation period or assimilation window, the background term (i.e., the last term in (6)) is less valuable since the situation of retrieval problem tends to be overdetermined. This tends to be the case for our 1D problem because of the relatively small degree of freedom. Also, for land surface variables, we usually lack proper background analyses a priori (e.g., usually no a priori soil moisture information). In this study, the measurement standard deviation is assumed to be 1 K and that for specific humidity is 1 g/kg. Apparently, the ratio of these two ’s plays a key role in the convergence of the retrieval procedure.
Suppose the vector form of the physical model (see (3a)–(3e)) is , where represents model parameters. Following the treatment in [22], one obtains the adjoint model asThe adjoint model, (7), is integrated backward in time to obtain ; the gradient of the objective function with respect to the initial conditions can then be computed asThe right-hand side of (7) determines how the information in and is assimilated. The second term on the right-hand side of (8) is trivial, since we have no proper prior estimate of the initial condition vector as argued previously.
We did not discretize (7) to obtain the adjoint system. Rather, we here follow a code-to-code approach [23, 24], which is especially suitable for sophisticated models with complicated boundary conditions and on/off switches. The adjoint model (ADM) constructed this way is also strictly consistent with the forward model. We also chose to hand-code the tangent linear model (TLM) and ADM because hand-coded adjoint is more humane and more efficient than machine produced one. More importantly, it helps the researcher to become more familiar with the operation of both the forward and backward systems, which comprises the prerequisite for good research. For example, in the forward scheme, the second-order Runge-Kutta method requires two calls of the subroutine tendency, which provides the right-hand side of the soil equations in [1]. In the adjoint coding, the synchronizing of these two half-time steps (for soil model integration) is necessary. We solved this dilemma by interim information storage.
Also, in the adjoint coding (not shown here for reasons of clarity), several continuous alternatives for stomatal resistance ’ parameterization are available to cope with the above-mentioned differential problems. Similar problem exists for the PBL parameterization for the bulk Richardson number. In calculating the bulk Richardson number, the convective contribution to the wind speed is omitted when the virtual potential temperature difference between the air and surface is less than 0.1 K. We adopted this criterion mainly because there is a singular point in the adjoint of the PBL component resulting from the parameterization of the convective component of wind speed. The gradient becomes infinite when the PBL experiences a transition from a stable to an unstable regime. We thus do not calculate the gradient for a temperature difference less than 0.1 K between the two temperatures used for the parameterization.
In formality, adjoint model is exactly a transform of the tangent linear model. The ADM code also is implemented strictly in the reverse order of that in TLM for the main driver routine and its subroutines. Using Fortran code as an example, in each statement, sensitivities of right-hand sided variables are expressed as the common left-hand sided variable. For example, if the TLM has a statement as , these become two corresponding statements in the ADM: and . The adjoint variables () accumulate the sensitivity information within the assimilation window and eventually become the gradients of the cost function with respect to the corresponding control variables (). Backward modeling (adjoint model construction), unlike the forward modeling (forecasting model and its tangent linear model), has no direct physical laws to guide the researcher; thus, finishing the adjoint code is merely half-way to finishing the adjoint system implementation. Thorough and objective testing of the components is indispensable. This involves the verification of the adjoint code and then the entire 4D-Var system.
Because ADM is based on TLM or, alternatively, the former is the matrix transform of the latter, the error sources of TLM are also the error sources of ADM. However, functions of TLM and ADM models are essentially different. TLM is used to approximate the evolution of a small perturbation on forward model. The time integration is forward. However, the backward integration of ADM is to unveil the sources of the perturbation. The adjoint variables corresponding to output of the forward model are actually inputs to the backward integration. Reflecting on the hand coding, for each statement, the perturbation of the right-hand terms is expressed in the form of perturbation on the left-hand side term.
To verify the correctness of the adjoint model, we use the inner product of the solutions to the direct and adjoint equations, a time invariant quantity [25] (Eqs. and in [25]). Consider that the magnitudes of different control variables are quite different. For example, soil moisture is less than 1.0, whereas soil temperatures are around 300. If all variables are perturbed simultaneously, the small quantities (considering only the numerical values, not the units) tend to be overshadowed by the larger quantities in the inner product. Therefore, we performed component-by-component verifications. Using double precision compiling, the relative difference is generally less than 10−8.
The function of ADM in a 4D-Var system is to estimate the gradient of the cost function with respect to the control variables. The correctness of the gradient value is verified against the Taylor expansion of cost function around specific control variables (p32 in [26]):where is a small real number, is a random perturbation vector which can be generated by using the Fortran library function, and is the transpose of .
For values of which are small but not very close to the machine zero, one should expect a value of approaching 1 linearly for a wide range of magnitudes of . We performed an experiment using coupled run of the land-atmosphere model and its adjoint with a cost function defined as in (6). The numerical results are shown in Figure 6. It is clearly seen that for values of between 101 and 10−4 unit values for are obtained. The curve deviates from unity for different reasons as and . The former is due to truncation error, while the latter is because the small perturbation assumption is violated. The correctness of the gradient of the cost function and the correctness of the vector product have therefore been verified. Interested readers can find a more sophisticated and efficient method in evaluating the TLMs using Taylor expansions [27].
4. Numerical Experiments
We use a variational method to determine the land surface model initial states for which the model forecast best fits screen level atmospheric observations of potential temperature and specific humidity within the assimilation window. This is achieved through minimizing a cost function (see (6)) defined as a quadratic difference of the screen level air temperature and specific humidity. is an implicit function of the control variable . contains the soil temperature ( and ), the soil water contents ( and ), and the canopy interception/dew . Because is insignificant in the selected periods, has only four effective dimensions. For OSSE (and actual measurements), the convergence criterion for the optimization procedure is when the cost function decreases by 3 () orders of magnitude or the iteration number exceeds 30 (50).
Using OSSE experiments, we discuss the retrieval scheme robustness to the magnitudes of initial guess errors and to magnitudes of Gaussian noise. In practical application, the shorter the assimilation window the better the performance. However, because of model error and measurement inaccuracy, sufficient information content is required for reliable retrieval. The optimal window length is apparently a function of sampling (data ingest) frequency. Some techniques for improving convergence of the optimization will also be addressed for OASIS data assimilation. With these issues in mind, we designed a set of numerical experiments (see the following).
Retrieval experiment design for coupled model system is as follows:
OSSE. We have a 3-hour nature couple run starting from 00 Z, 6 August 2000; synthetic (simulated) observations of atmospheric potential temperature (PT) and specific humidity () were sampled every 5 minutes and used by the retrieval experiments; = (295.8 K, 298.5 K, 0.253 m3m−3, 0.278 m3m−3, 0.0 m2m−2):(a)Initial guess errors exist in all control variables: = (300.8 K, 303.5 K, 0.303 m3m−3, 0.328 m3m−3, 0.0 m2m−2) or = (5 K, 5 K, 0.05 m3m−3, 0.05 m3m−3, 0.0 m2m−2) (OSSEa).(b) in OSSEa was halved in addition to several other representative experiments (OSSEb).(c)There is scheme resistance to Gaussian noise in PT and (OSSEc):(1)Gaussian noise series of zero mean and different standard deviations added to obs.(2).(3).(4)Both and (guarantee that noise series added on and are not correlated).
OASIS. OASIS observations of PT and available every 5 minutes are used by the retrieval experiments:(a)Control run: initial guess control variable the same as OSSE runs (OASISa).(b)Effects of data availability (OASISb):(1)data frequency: 5 min, 30 min, 1 hr, 3 hr, and 6 hr;(2)data availability period: running over the daily cycle with hourly interval and with assimilation window length from 3 to 24 hrs.(c)Identifying the optimal assimilation window length (OASISc).
4.1. Observing Simulation System Experiment (OSSE) Retrieval
We simulate observations every minute of temperature and specific humidity at 2 m. The retrieval scheme may use the data less frequently. For the first retrieval experiment (OSSEa), we assume initial guess errors in all control variables (i.e., , all significant compared to typical existing model errors for such quantities. The reference run is a coupled 3-hour nature run starting from 12 Z on 6 August 2000 with true initial states: , and . Synthetic observations of atmospheric potential temperature and specific humidity are sampled every five minutes and used by the retrieval experiments.
The first set of experiments is performed assuming perfect observations. The retrieval procedure is shown in Figure 7 within the assimilation window. We see that the retrieval method is able to realize a good analysis, that is, to approach closely the real initial states in just seven steps. The retrieved initial values are . Compared with the true values, the differences are insignificant , , (Figure 7(a)). As a result (of the successful adjustments to the initial guess errors on the states) the cost function decreased significantly to about 2.54 × 10−4 at step seven (Figure 7(b)). The trajectory fitness within the assimilation window is overall very satisfactory for all four state variables (Figure 7(c)) and surface fluxes (Figure 7(d)). The initial hot moist guess signifies a drastic overestimation to LE but underestimation to . The retrieved curves are indiscernible to the synthetic true trajectories.
(a)
(b)
(c)
(d)
Using the retrieved initial states, we made extended forecasting of 24 hours (Figure 8). The improvements of using retrieved over initial guessed states are significant, especially the case for daytime heating period (12–24 Z), as measured by the statistics as summarized in Tables 2 and 3. MAE means maximum difference of two time series; MBE means mean bias error; and RMS means root mean biased error. For example, the RMS errors for screen level air temperature, specific humidity, LE, and are reduced from 3.39 K to 0.15 K, from 6.22 g/kg to 0.05 g/kg, from 17.42 W m−2 to 0.62 W m−2, and from 13.2 W m−2 to 0.52 W m−2, respectively. The improvements resulting from using retrieved over initial guess state variables are also shown in the trajectory fitness of the state variables themselves (columns 4–7 in Tables 2 and 3).
4.1.1. Robustness to Initial Guess Magnitudes
In addition to the case discussed above, we also tried using the initial guess control variables halved (the second experiment listed in Table 4). The retrieval is a very successful one in that the retrieved states are close enough to the true states (). Retrieval from using further reduced initial guess state variables (e.g., the third experiment listed in Table 4) is always successful in accordance with our tangent linearization assumption. The fourth experiment listed in Table 4 is a rather free guess of 300 K soil temperatures and 35% wetness on soil moisture. The corresponding retrieval is a successful one. Actually, the scheme is very robust for initial guess on soil temperature (270–320 K). To investigate the scheme robustness to soil moisture initial guess errors, we performed two retrieval experiments with the same soil temperature initial guess of 300 K (the last two as listed in Table 4). For the very dry initial guess case, the adjustment to deep soil moisture is insignificant, whereas for the wet initial guess, the initial guess error in surface soil moisture cannot be effectively removed. This can be explained through analyzing the structure of the cost function (note: its production is rather time consuming), especially the cross section.
Around and below wilting point, the isoligns tend to be vertical and the cost function is no longer sensitive to initial deep soil moisture value. Above saturation point, the isoligns tend to orient horizontally and the cost function becomes insensitive to initial value of surface soil moisture values. Generally speaking, for this relatively dry period, the separation of sensitivity between and is not so strong that it is easy to get both quantities successfully retrieved.
4.1.2. Robustness to Gaussian Noise
Another important issue is the noise level that the scheme can tolerate. To this end, the observations created with the coupled 1D land surface atmospheric model are corrupted with realizations of a Gaussian noise with zero mean and different standard deviation magnitudes on and . Using the reference retrieval experiment configuration (OSSEa), we performed retrieval experiments with Gaussian noise of different magnitudes (0.1, 0.5, and 2.0 K for temperature and 1, 2, and 5 g/kg for specific humidity) added to the air temperature and mixing ratio time series within the assimilation window, respectively, and jointly. We compared retrieval experiments with noise added to both screen level parameters and those without noise corruption. Generally, the higher the noise level, the worse the retrieved initial states and consequently trajectory fitness (as measured by the cost function). Also, the retrieval of deep soil moisture and soil temperature seems to be affected more severely by the corrupted air temperature observations, whereas the bad effects from corrupted specific humidity observations are felt by all control variables. For Gaussian noise added to the 2 m air temperature, when the standard deviation is smaller than 4 K, convergence is not apparently affected. For Gaussian noise added to specific humidity, standard deviation smaller than 3 g/kg has no apparent effects on convergence.
4.2. Retrieval Experiments with Real Data
The first real data assimilation experiment (see Figure 9) uses exactly the same initial guess as OSSEa but instead assimilate the OASIS observed screen level air potential temperature and specific humidity. Similar to the corresponding OSSE experiment (Figure 7), all control variables get significantly improved (Figure 9(a)). Primarily, due to the model errors, the decrease of the cost function is not as apparent as the corresponding OSSE. Upon convergence, there is still ~36% error as measured by the cost function (Figure 9(b)). Meanwhile the trajectories of the prognostic variables are significantly updated (Figure 9(c)) within the assimilation window due to the accurate retrieval of the control variables (. Consequently, the improvement in surface fluxes is significant (Figure 9(d)). We further examine the improvements for the one-day long extra forecasting period that follows the assimilation window period (Figure 10). The improvements are all apparent and significant, especially for the daytime heating period. The RMS error for is reduced from 13.2 W m−2 to <5 W m−2. For LE, RMS error is reduced from 17.3 W m−2 to merely 4.7 W m−2. More detailed statistics are given in Tables 5 and 6.
(a)
(b)
(c)
(d)
Unlike the corresponding OSSE, the retrieval of is good only on the daily or longer time scale basis. Although the mean value and the amplitude are correctly retrieved, there is still significant phase error associated with the retrieved time series. As stated above, forward model needs to include more sophisticated vegetation processes to improve the simulation of the daily cycle.
4.2.1. Influence of the Assimilation Window Length
For OSSE experiments without noise corruption to the simulated observations, there is no apparent benefit in using longer assimilation window length (>6 hours for 5 min ingesting frequency). Retrieval experiments with real OASIS observations are performed for six assimilation window lengths: 1, 2, 6, 12, 18, and 24 hours. The performance of the retrieval system, as measured by the closeness of the retrieved initial states to true states (), can be seen from Table 7. For assimilation window less than 12 hours, increasing of assimilation window is beneficial. Choosing assimilation window longer than 3 hours can steadily reduce the initial guessed error by more than one order of magnitude (considered significant reduction for real data assimilation). Further increasing the assimilation window length from 6 to 12 hours is less effective in further reducing the initial guess error. More interestingly, using too long assimilation window (e.g., 18 hours) sacrifices the retrieval. This degrading trend is consistent from 12 hr–36 hr assimilation window length and becomes not practical beyond 36-hour assimilation window length. This signifies contesting between model error and the information content.
Another severe hindrance of using too long assimilation window is the computational cost associated with both forward model integration and the backward adjoint integration. For the sake of efficiency, we suggest assimilation window of about six hours for this relatively moist period and vegetated surfaces. From Table 7, for this data ingesting frequency (5 minutes), using 6-hour assimilation seems a reasonable choice. The retrieved temperature values can be within 1 degree of the true (0.75 and 0.95 K, resp.) and the retrieved soil moisture is even close (0.005 and 0.004 m3 m−3, resp.).
For the shared forecasting period of 12 Z 7 August-12 Z 8 August, we compare three assimilation window lengths: 1 hr, 6 hr, and 24 hr. The 6 hr assimilation window gives the most desirable retrieval results in this case. The RMS errors are 0.42 K for air temperature time series as compared with the OASIS observations, 6.35 × 10−4 kg kg−1 for specific humidity, 0.83 K for skin temperature, 0.23 K for deep soil temperature, 7.89 × 10−3 m3 m−3 for superficial soil moisture, 1.08 × 10−2 m3 m−3 for deep soil moisture, 11.7 W m−2 for sensible heat flux, and 15.6 W m−2 for latent heat flux, respectively. The RMS errors using too short (1-hour) assimilation window are much worse, 3.29 K for air temperature, 3.09 × 10−3 kg kg−1 for specific humidity, 5.09 K for skin temperature, 2.08 K for deep soil temperature, 3.81 × 10−2 m3 m−3 for superficial soil moisture, 3.42 × 10−2 m3 m−3 for deep soil moisture, 61.1 W m−2 for sensible heat flux, and 84.2 W m−2 for latent heat flux, and so are those using too large assimilation window length (24-hour), 0.61 K for air temperature, 9.37 × 10−4 kg kg−1 for specific humidity, 1.05 K for skin temperature, 0.47 K for deep soil temperature, 1.16 × 10−2 m3 m−3 for superficial soil moisture, 2.37 × 10−2 m3 m−3 for deep soil moisture, 19.1 W m−2 for sensible heat flux, and 24.3 W m−2 for latent heat flux. Other aspects of the comparison can be seen from Figure 11. In Figure 11, we plotted the time series of all the control variables, the latent and sensible heat fluxes, and the screen level air temperature and humidity resulting from using the final retrieved initial land surface conditions. Measured by all three statistical indices (MAE, MBE, and RMS), the 6-hour assimilation window length performs the best. For example, the mean bias error for LE is 5.3 W m−2, significantly less than 62.7 W m−2 from using a 1-hour assimilation window and 20 W m−2 from using a 24-hour assimilation window. The maximum difference is only 45 W m−2 for LE, whereas it can be as large as 190 W m−2 when using a 1-hour assimilation window length.
Using assimilation window lengths from 4 to 24 hours and initial guess control variables as in the control run, we addressed another related issue about influences of the assimilation frequency. For each assimilation window length, we tested data ingestion frequency of 5–60 minutes with increment of 5 minutes. A subset of this experiment (assimilating the observations every 10, 15, 30, and 60 minutes) is listed in Table 8. For each combination of assimilation window length and data ingesting frequency, we provided the final retrieved surface temperature and moisture values as indices for evaluating the performance of the retrieval.
As expected, for every assimilation window length, less frequent data ingesting by the adjoint integration is much more efficient. Generally speaking, the longer the assimilation window is, the more deserving to perform less frequent assimilation. For example, 18-hour assimilation window can tolerate assimilating the data every 60 minutes, whereas 4-hour assimilation can at best assimilate data grid every 10 minutes. When assimilation frequency becomes too low (<every 10 minutes in this case), the retrieved states will deviate from the true states significantly and even fail to converge. The apparent benefit for performing assimilation at lesser frequency is for efficiency of the scheme. Updating the forcing every minute makes the scheme less practicable for real application. For example, for some initial guesses, it takes 2 hours for a single column run of the retrieval scheme for a 24-hour assimilation window and ingesting the observations every 5 minutes.
For a given assimilation window length, it seems that the more frequent data ingesting yields better retrieval. However, as assimilation becomes less frequent and the assimilation is long, the uncertainty increases. For example, for 18-hour window length, using hourly data ( K) performs better than using 30-minute data ( K).
As mentioned above, as data ingesting becomes less frequent, the optimal assimilation window length tends to increase. For example, 15-minute data supports an optimal window length of around 18 hours. There is no edge for using window lengths between 18 and 24 hours since all retrieval experiments arrive close to K, K, m3 m−3, and m3 m−3.
An interesting aspect we would like to mention is that the degrading effects on retrieval from less frequent data availability are quite different for different components of the control variables. For shorter assimilation windows (e.g., 4-hour and 6-hour in Table 8), the retrieval of is deteriorating fast (298.65 K using 5-minute data and 300.79 K when using 10-minute data for 4-hour window length). However, for longer assimilation window lengths (e.g., 24-hour), the data frequency has no major effects on the retrieval of surface soil moisture value. This can be explained by analyzing the sensitivities of screen level air potential temperature and humidity time series to variations in initial and (Figure 12). In Figures 12(a) and 12(b), initial value of is varied down to 0.1526 m3 m−3 and up to 0.3526 m3 m−3, with increment of 0.05 m3 m−3. Figure 12(a) indicates clearly that the responses from screen level potential temperature are significant at all phases of a daily cycle, whereas the response from specific humidity to initial value of is more salient during nighttime after 18 Z and is insignificant at the beginning 6 hours. To see the sensitivity of screen level atmosphere to initial value of , we provided Figures 12(c) and 12(d), where initial value of is varied down to 285.79 K and up to 305.79 K, with increment of 5 K. In contrast to the sensitivities to , all sensitivities are contrasted in the first couple hours for variations in initial value of . This explains why decreasing the assimilation frequency sacrifices the retrieval of initial , especially when window length is short (for 4-hour window length, 298.65 K versus 300.79 K in Table 8). This also partly explains why the retrieval of initial is always pretty accurate for longer assimilation windows.
(a)
(b)
(c)
(d)
We also find that there is little benefit in using assimilation window longer than 24 hours. Except the apparent reason of model inaccuracy, nonlinearity may also play a role in scarifying the usefulness of the adjoint technique, which is based on tangent linearization as mentioned above.
4.2.2. Sampling Strategy
The general setting is identical to the OASIS control run. For the following set of retrieval experiments, we, however, assume data availability only at certain periods of a daily cycle. We tested assimilation window length from 1 hour to 12 hours, and for each window length we tested centering it at different times of a day. Results from applying different window lengths unanimously point to the daytime period being more informative for the state retrieval purpose. This finding is in accordance with [2]. We also find that using a window length longer than 4 hours and centered at local noon can be as informative as using data within all the daily cycle for retrieval of the land surface states. Certainly, this critical window length of 4 hours is connected with the relatively high data ingesting frequency (5 minutes each). We thus tested less frequent data ingesting and enlarged assimilation window length (and the overdetermination requirement is always satisfied) and find that this above assertation still holds qualitatively; that is, there exists a window length centered at local noon significantly shorter than 24 hours and is nearly equally informative as the whole daily cycle period.
We now understand better why the retrieval significantly improved when using assimilation window longer than three hours for our experiments testing effects from assimilation window length. The assimilation windows all start at 12 Z, 6 August 2000. It is a heating up stage of the daily cycle. Except that longer assimilation window length contains more observations, using assimilation of 1 hour is also closer to the nighttime period than using, say, 4-hour assimilation window.
5. Summary and Conclusions
Modern numerical weather prediction requires the initial conditions of the land surface/soil-vegetation model. The lack of routinely available measurements of land surface variables motivated the development of schemes that infer land surface information from routinely measured variables.
Based on 4D-Var strategy, a retrieval system for the initial state of a soil-vegetation model is developed which assimilates operationally available screen level atmospheric measurements. Coupled soil-vegetation and the atmospheric models and their adjoint have to be constructed. To avoid the need for developing a full 4D-Var system of a 3D coupled system, a 1D column simplifying assumption is made. Radiation and vertical mixing are the main processes involved. The vertical mixing is described using the nonlocal MRF PBL scheme, and soil-vegetation processes are described by a two-layer force-restore model implemented in ARPS.
The sensitivities of the measured variables to the retrieved (control) variables are first studied via forward simulations. Such simulations are verified against OASIS measurements. An important correction is made to the deep layer soil temperature prediction equation to enable the successful retrieval of the deep layer temperature and improve that of others.
Unlike previous studies [2, 8], both OSSE and real data retrieval experiments were performed, assuming two different types (skin temperature or screen level atmospheric variables) of measurements, for a dry period during year 2000 at the Norman OK OASIS site. The retrieval system is able to recover the correct values of soil state variables in most cases. The robustness of the systems is assessed by testing different magnitudes of initial guess errors and different measurement errors. The impact of the assimilation window length and data availability on the retrieval products is also examined.
Compared with [28] and several other oversimplified treatments of land surface processes [2, 29], our explicit treatment of vegetation processes is critical to the successful real data assimilation for deep soil moisture content for partially vegetated surfaces. In contrast, Xu and Zhou [29] discussed a linear regression method for retrieving bulk soil moisture contents from soil temperature measurements, based solely on the soil heat capacity’s dependency on soil moisture contents. Results from [2] are also limited to the bare ground case.
Our retrieval system does not perform effectively for very wet circumstances. Ren [1] performed a series of retrieval experiments for such conditions. The most important finding by Ren [1] was that the measurement (cost function) was shown to be insensitive to the value of deep soil moisture for very wet conditions. The marked sensitivity contrast between and prevents all control variables from being simultaneously retrieved successfully without proper preconditioning.
In this study, by restricting the application of our data assimilation system to only weather situations when impact of soil moisture on near surface atmospheric conditions is dominant, for example, to situations with clear sky, low wind speed, and strong radiative input, it is certain that the analyzed soil moisture/temperature corrections are not caused by an attempt to compensate for model deficiencies in parameterization. For more sophisticated systems, the analysis will be more difficult to obtain. Significant efforts will be put into that direction in the future.
In reality, the atmospheric measurements are available less frequently (usually 6 hours for SYNOPs) than required by our retrieval scheme. Further improvement to the land surface scheme is required to reduce the phase error in the simulation of the specific humidity daily cycle. Through forward model improvements, we believe the robustness of the retrieval scheme when supplied with real observations will be improved and can work more effectively with less frequent observations.
The adjoint of a much more sophisticated land surface scheme, the Community Land Model (CLM 3.5), is now available [30]. The present study, performed a decade ago, still has significant value in guiding the operational use of adjoint-based variational schemes such as the variational data assimilation system of CLM, specifically in the specification of the assimilation window (its temporal location and length) and the ingest frequency for observational data. Although computationally very efficient, the adjoint-based variational data assimilation procedure involves coding of the adjoint model, which must strictly correspond to the forward model. Subsequent changes made to forward model require essential modification to the adjoint model. This is a technical disadvantage and explains the prevalence of other forms of data assimilation and optimal parameterization, especially the Ensemble Kalman Filter [31–37]. For these recent studies, our results on information redundancy and the observational error covariance provide useful guidance for such operational applications.
Appendix
For any quantity (can be , , , or ), the nonlocal approach for vertical mixing, based on [18], can be expressed asThis equation is a form of closure of the diffusion (first term) term in (1a)–(1d). Since local mixing uses to estimate vertical turbulent fluxes, under conditions that the mean quantity is the same for the two selected level, zero flux will result. Chances are actual bulk fluxes are not zero (may flip signs between the two levels). Equation (A.1) is proposed also to deal with this shortcoming with local schemes. This also suggested that parameter should be a bulk property of the total convective boundary layer not that of a quantity which represents a subdomain of it. As in [16, 18, 38, 39], the diffusivity coefficient for momentum is formulated aswhere is Karman constant, is a characteristic turbulent velocity scale for scalar transport, and is the depth of the PBL. The mixed layer velocity scale is represented asFollowing formulation of [40] (Eqs. ), here is evaluated at the top of the surface layer , with the stability criteria expressed as the covariance of virtual potential temperature perturbation and vertical velocity evaluated at the surface. For example, signifies unstable condition. The PBL depth in (A.2) is iteratively obtained throughwhere is critical Richardson number (set to 0.3 in this study); is the virtual potential temperature at the lowest model level and should be taken consistent with in the following definition of bulk Richardson number:Here is the appropriate temperature near the surface (about 10 m, or near the top of the surface layer as in [38]) defined asThe right-hand side of (A.6) is the scaled virtual temperature excess near the surface, usually bounded (3 K in [16]) in case surface wind is very weak. Here is a coefficient of proportionality (taken as 7.8 in [16]), which is used in describing the nonlocal enhancement effects from large-scale eddies to total flux. ConsiderSimilarly, the dimensionless vertical temperature gradient given by [40] isThe eddy diffusivity for temperature () and moisture () is estimated from the relationship of the turbulent Prandtl number:Following [16], for the free atmosphere, a local- approach [41] is utilized: , where the mixing length ; here is the asymptotic length scale, and is the stability function. The stability function is represented in terms of the local gradient Richardson number at a given level. Here, is the virtual potential temperature. Computed is usually numerically bounded to prevent unrealistically unstable regimes (−100 in [16]).
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work was supported by DOT-FAA Grant no. NA17RJ1227-01 and NSF Grants nos. ATM9909007 and ATM0129892. The authors thank Dr. J. Brotzge for making available the OASIS data set and for very helpful discussions on the related issues. The lead author’s Ph.D. thesis under Professor M. Xue forms the basis of this research. Professor Mervyn Lynch (Curtin University) proofread the entire paper.