Evaluating the Technique of Variational Assimilation in a Semi-Arid Environment
by
Karen I. Mohr
Class project for CE 394K GIS in Water Resources
c. Preparing the Soil and Vegetation Database using GIS
Soil moisture has a significant influence on the continental surface energy budget. One of its important roles is the partitioning of energy fluxes into latent and sensible heating. This partitioning is believed to exert a major influence on the development of convective clouds and precipitation. Unfortunately, soil moisture observations are scanty to non-existent, often forcing weather and climate modelers to assume climatological, even arbitrary, values for soil moisture in their experiments. This project initiates a study to assess the usefulness of a technique proposed by Mahfouf (1991) for calculating soil moisture values using a land-surface model driven by atmospheric observations. Study experiments will be conducted in the Little Washita, OK watershed, a small watershed in a semi-arid environment..
In the initial stages of the study, GIS was an excellent tool for organizing and evaluating spatial data required for assimilation into land-surface models. Three major tasks were performed in GIS. Soil and vegetation data were easily linked to the grid cells of the model domain for assimilation into the land-surface model. Index files specifying the properties of soil types and land cover in the model domain were developed from GIS attribute tables. Required atmospheric data from widely, irregularly spaced weather stations were interpolated to the model grid.
Because of unforseen delays, results of modeling experiments cannot be presented here, and this project report treats the study as a work in progress. In the future when experimental results are available, it appears GIS will also be useful in evaluating these results spatially.
A number of studies have demonstrated the sensitivity of continental weather and climate dynamics to soil moisture. Sellers et al. (1989) identify three main land-surface properties governing land/atmosphere interaction: surface albedo (radiative transfer), surface roughness (momentum transfer), and surface soil moisture (latent and sensible heat transfer). Soil moisture is directly involved in determining the magnitude of surface heat transfer and its partitioning into latent and sensible heating. It is also indirectly involved in the former two properties.
Unlike atmospheric temperature, pressure, wind speed, and humidity, soil moisture is not a commonly observed geophysical variable. Soil moisture observations are scanty to non-existent except for limited time periods in limited areas during major field experiments. Maps of soil moisture are produced either from remote sensing, climatology, or model-generated. For large-scale climate models, the natural variability of soil moisture can be reduced to a few harmonics, obviating the need for precise spatial representation of soil moisture at less than monthly time scales. This is not a safe assumption with 100 km-scale weather models operating at hourly time scales. Hence, Mahfouf (1991) describes a technique, called "variational assimilation", to generate variable soil moisture fields for land-surface models running at short time scales.
The "variational" part of the title refers to the use (i.e. assimilation) of time-varying atmospheric data. Variational assimilation uses common surface meteorological observations such as temperature and relative humidity and climatological soil moisture values as input to a land-surface model . The land surface model is run until equilibrium is reached between evaporation and precipitation, typically 48 hours at the mesoscale. The soil moisture values evolve as water balance is achieved, such that the final soil moisture patterns reflect the final atmospheric conditions. Comparing model results to soil moisture data from a couple of case days in the HAPEX-MOBILHY field program, Mahfouf (1991) suggests that realistic soil moisture patterns can be achieved using this technique. While variational assimilation is not preferable to actual data, it may be an improvement over a priori assumptions. These model-generated soil moisture patterns can then be assimilated into weather models to study land/atmosphere interaction in the formation of clouds and precipitation.
This study is not GIS-based but looks at where it may be GIS-related. The goal of this study is to see how well the technique of variational assimilation can derive realistic soil moisture patterns in a small watershed in a semi-arid environment during a dry down. It is the objective of this class project, which represents the initial stage of the study, to see where and how GIS may be integrated into the experimental process. Since large amounts of spatially distributed data are required by the land-surface model chosen for this project, it seems obvious that GIS could be useful. Where and how much are the questions to be answered.
NB: A related study by Seann Reed: A Soil-water Balance and Continuous Streamflow Simulation Model that Uses Spatial Data from a Geographic Information System (GIS)
a. Study area
Some large field experiments collect soil moisture observations. Since in situ soil moisture observations require a significant amount of human labor, temporal and spatial resolution are limited. The Southern Great Plains 1997 Hydrology Experiment (SGP) is unique in its collecting of soil moisture observations in a mesoscale domain over several weeks. This time period (mid-June to mid-July) encompassed two major precipitation events in which the evolution of soil moisture during a dry-down could be studied. Since evaluating the technique of variational assimilation requires ground truth soil moisture observations for verification purposes, this study looks at the Little Washita, OK watershed during SGP, 9-16 July 1997.
The Little Washita (LW) watershed in southwestern Oklahoma, straddling Grady, Caddo, and Comanche counties. The LW is part of the Washita River (HUC 111303) system, in the Upper Washita section (HUC 11130302). The Upper Washita is outlined in green, the LW in yellow. The map is a combination of the counties coverage developed in Exercise 3, and the HUC and RF1 files from Exercise 4. The coverages were re-projected in Arc/Info to a common Albers EA projection.
The LW summer climate is semi-arid, characterized by infrequent showers followed by days of dry, hot weather. Precipitation in this environment tends to be highly variable in both time and space, making it a good testing ground for the technique of variational assimilation. The time period (9-16 July) was chosen because a significant precipitation event occurred on 11 July. This time period covers both pre-event conditions and the post-event dry down.
The model chosen for this study is the Parameterization for Land-Atmosphere-Cloud Exchange (PLACE) developed by Peter Wetzel and Aaron Boone of NASA-Goddard Space Flight Center, MD. The PLACE model is a detailed land-surface process model of the partly cloudy atmospheric boundary layer and underlying heterogeneous land surfaces (Wetzel and Boone 1995). The vertical sub-surface column maintains five reservoirs for liquid water from the surface interception store down to bedrock. There are seven interactive layers for ground heat storage and exchange. PLACE consists of linked process models (e.g. net radiation, turbulent sensible and latent heat fluxes, ground heat storage) to predict surface temperature and boundary layer depth and growth. Turbulent eddy formation, vegetation responses, and soil water fluxes evolve from their initial state in response to the surface energy balance, closing the energy feedback loop. PLACE code is written in FORTRAN-77.
Click here to see the schematic of PLACE more clearly
PLACE requires a large number of data inputs. Click here to see the data dictionary. The model reads four index files to specify soil classes and their properties, vegetation classes and their properties, model constants such as the critical Reynolds number, and initial soil moisture and temperature. The data dictionary hints at the difficulty of assembling a database for this study. Although PLACE is intended for mesoscale, short-term applications, PLACE is a scale-invariant model, allowing the user to define a grid as fine or as coarse over as large a domain as needed. Since the data grids obtained for this study were 1 km x 1 km resolution, the model grid cells were also specified to be this resolution. Because of the computational overhead involved in such a fine mesh--PLACE has interactive grid cells--the size of the domain was limited to 50 x 50 1-square kilometer grid cells. During experiments, each of the 2500 grid cells is assigned a soil and land cover class from which PLACE calculates the energy fluxes and water storage.
c. Preparing the soil and vegetation database using GIS
The SGP experiment collected in situ gravimetric soil moisture and bulk density measurements at 21 fields around the LW watershed. A network of 42 weather stations, a micronet, were installed by the USDA Agricultural Research Service in 5- to 10-km intervals around the watershed, and there were several already installed weather stations belonging to the Oklahoma state mesonet network. These locations were transformed into point coverages and are plotted on the DEM below.
The atmospheric data required to run PLACE were derived from the mesonet stations. Although the micronet and mesonet stations both collect atmospheric data in 5-minute increments, the micronet stations do not measure wind speed or pressure, two atmospheric variables required by PLACE. Thus, data from the micronet stations have not been used at this point in the study. The LW fields in the watershed were joined to tables containing bulk density and gravimetric soil moisture measurements for the study time period. These observations could also be plotted on the map if desired. Click here to see the soil moisture data table. Eventually, the soil moisture observations will be used in verifying model generated soil moisture fields.
The DEM above is a 100-m resolution DEM and is only used here for display purposes. The grid (e.g. DEM, land use/land cover) and polygon (e.g. USGS STATSGO) data used to develop the model database were, as stated earlier, 1 km x 1 km resolution. The source of these data is the Earth System Science Center at Penn State University. The Earth System Science Center has assembled spatial data coverages in Arc/Info export format for the SGP region. The soils data are from the USGS STATSGO database. Below left is a map of the model domain and the LW watershed superimposed on the STATSGO polygons in the LW region. The red gridlines outlining the model domain were created in Arc/Info using the GENERATE function with the FISHNET option. The fishnet of polygons is created from specifying the x-y origin and the number of rows and columns from the origin.
(Above left) Model domain and STATSGO polygons. (Above right) Model domain and land cover/land use grid. Major land cover types: purple = woody savannah, orange = grasslands, blue = crops, green = mixed crops and natural vegetation, pink = evergreen needleleaf forest.
Within the model domain are 15 different MUID and 8 different land cover types. The land cover/land use grid shown is for the International Geosphere Biosphere Programme (IGBP) global land cover categories. These are consistent with the USGS Anderson classification codes, but the IGBP codes are more specific with vegetation types, making it easier to assign the correct vegetation properties (e.g. stomatal functioning temperatures) to each grid cell in the model domain. The union of the model domain polygons and the STATSGO polygons in Arc/Info created a table of model grid cells with their corresponding STATSGO MUID. Similarly, the land cover classes can be unioned to model grid cells. These tables were imported into Excel, assigned an index (1-15) for soil and (1-8) for land cover, converted to a text file, and then this initialization file can be assimilated into PLACE to assign soil and land cover classes to each model grid cell.
Index files read by PLACE assign values to the variables in the data dictionary for each soil and land cover class. Among the Earth System Science Center's SGP data files are INFO tables with vegetation types and percent cover referenced by MUID and component. These could be joined to MUID component and layer tables making a complete description of soil and land cover characteristics for the MUIDs in the model domain. Click here for an example. Several assumptions/simplifications were made in creating the index files from the spatial data tables. Soil properties, particularly sand/silt/clay fractions, are averages of the properties of the largest components in an MUID. The Excel workbook in the image above is where the INFO tables were exported to perform these calculations. The averaging assumes that an MUID's major components are present in the LW and in all the grid cells assigned that MUID in similar proportions. A similar assumption is required with plant cover. Vegetation properties (e.g. stomatal resistance) were drawn from general reference books (e.g. Montieth and Unsworth 1990) and are based on class rather than species. They are also averaged for the major plant types in a land cover category. Since the INFO tables list actual species, further research in the future may be able to improve the representativeness of these properties. Below is an example of the vegetation index file. The numbers in the first column refer to the variables listed in of the data dictionary.Section 3
Class 7 is surface water. PLACE allows grid cells to be defined as water bodies, and it adjusts the calculation of evaporation rates accordingly.
Two of the soil properties concern topographic changes, delta topography, and delta horizontal. These properties can be thought of as a sine wave, with delta horizontal referring to the distance between wave crests and delta topography referring to the amplitude of the wave. Overlaying the STATSGO polygons on top of the DEM (below) suggested that each MUID in the model domain could be assigned a delta topography in tens of meters and a delta horizontal in thousands of meters.
DEM for model domain overlaid by STATSGO polygons
Both of these properties are crude approximations for each MUID in the model domain. Since PLACE has simple parameterizations for horizontal flow, further precision was deemed unecessary at this point in the study.
d. Preparing the atmospheric forcing data using GIS
Section 4 of the data dictionary lists atmospheric forcing variables PLACE requires to calculate surface heat and water fluxes. To obtain these values PLACE can be coupled to an atmospheric model or can be run in stand-alone mode, in which the user specifies a file of atmospheric data. An example is shown below. Column names are, left to right, shortwave radiation, longwave radiation, temperature, wind speed, precipitation, station pressure, and specific humidity.
At each model time step, PLACE reads a line of atmospheric data corresponding to each grid cell. For 2500 grid cells, 2500 lines are read at each time step, in this case, 5 minutes (300 seconds). The time step here is the same as the measurement interval for the OK mesonet stations. The number of lines in the atmospheric data file for a simulation time of 1 day is (1 day * 86400 s/day * #grid cells)/time step in seconds. For 2500 grid cells, this corresponds to 720,000 lines of data.
The raw mesonet data required extensive processing to put it in the form shown above. Several rounds of data processing took place. The first involved extracting the sites of interest from the statewide network data files and putting the desired sites and desired variables into a common file by date. The next step was to calculate longwave radiation and specific humidity from observed values of temperature, pressure, and relative humidity. The last and biggest step involved taking irregularly spaced mesonet sites and interpolating the values observed at the mesonet sites to the model grid. PLACE requires specific humidity and both downward long- and shortwave radiation, but the mesonet sites measure only relative humidity and shortwave radiation. From the air temperature and station pressure, the specific humidity and longwave radiation can be calculated (formulas from Rogers and Yau 1989 and Dingman 1996):
Two simplifications were made in calculating the longwave radiation. First, there was no forest canopy. Since there are few forested areas in the model domain, for most of the domain this is a safe assumption. Second, a simple algorithm was developed to allow the fractional cloud cover to vary with time of day and with the presence of precipitation as a linear function. Climatologically, diurnal cloud cover in southwestern OK is a maximum at sunrise and a minimum at sunset. The calculation of C assumes 0.6 cloud cover at sunrise, decreasing in hourly increments to 0.2 cloud cover at sunset. The algorithm adds 0.6 to the normal cloud cover (to a maximum of 1.0) for six hours after the onset of recorded rain. Variations in downward longwave radiation from cloud cover and forest are in the tens of W/m^2, but it is not at the present time possible to see the effect error in the calculation of longwave radiation has on the calculation of latent and sensible heat fluxes by PLACE.
There are three mesonet stations inside the LW, and one just outside of it. Since the model domain extends for some distance north and south of the LW, other mesonet stations were included in the atmospheric data. Below left is a map showing all of the stations used and their positions relative to the model domain.
(Above left) map of mesonet stations relative to model domain and Washita River system. (Above right) the mesonet station in Ninnekah, OK in the northeast corner of the LW watershed.
A total of nine irregularly spaced mesonet stations were used to interpolate values for the atmospheric variables in the model grid cells. They are (inside the grid, clockwise from top right) Chickasha, Ninnekah, Acme, Apache, and (outside the grid, clockwise from north) Minco, Washington, Ketchum Ranch, Medicine Park, Ft Cobb. Because of ArcView's limited number of interpolating functions and the large amount of data involved, the interpolation was performed in IDL. IDL is a sophisticated graphics programming environment for analysis and visualization of many types of scientific data. Although IDL is not marketed specifically as a GIS tool, spatial data analysis applications can be written in it using any of the 14 different gridding and interpolation schemes. For this study, the atmospheric data from the mesonet were first fitted to a regular grid by Delaunay triangulation and then both spline-type and linear interpolation routines were tested. Spline interpolation was immediately eliminated as an option because extrapolation is inherent in the routine, and it produced negative values of precipitation. A linear interpolation scheme (Cressman weighting) with no extrapolation was chosen. Below is an example of a surface of precipitation data generated for one 5 minute time step:
This image is from 11 July 2355Z. A series of heavy thunderstorms thoroughly soaked the LW watershed earlier in the day. Note that all of the precipitation values (z-axis) are non-negative.
The assembling of the database was a time consuming process to gather and modify the necessary coverages. The atmospheric data required extensive processing on a powerful Silicon Graphics workstation. Several programs in C and IDL were written to accomplish the atmospheric data processing. All these tasks were completed in a timely fashion, but did not allow much time for actually conducting experiments on PLACE.
PLACE is composed of eleven subroutines linked to a driver module. The code distributed to the public by NASA-Goddard Space Flight Center includes sample files to allow the user to compile and test the code. Unfortunately, although the code comes with documentation on the physics and development of PLACE, it has almost no documentation on how to set up and run the model. The public code is hardwired to run the sample files which specify a domain of a single grid cell. In the limited time left after data processing, extensive modification was required to compile PLACE to assimilate this study's database. Since modifications were completed, additional problems, principally bugs in the subroutine calculating ground heat flux, have been identified. Since the ground heat flux subroutine is the most complex of the subroutines, debugging it will prevent this study from obtaining results in the class project time frame.
This project initiated a study to assess the usefulness of a technique, "variational assimilation", proposed by Mahfouf (1991) for calculating soil moisture values using a land-surface model driven by atmospheric observations. "Variational" refers to the use (i.e. assimilation) of time-varying atmospheric data.to generate variable soil moisture fields for land-surface models running at short time scales. This work was driven by the need to have realistic soil moisture patterns for studying land/atmosphere interaction. The goal of this study is to see how well variational assimilation can perform in a small watershed in a semi-arid environment during a dry down.
The land-surface model chosen for this study is the Parameterization for Land-Atmosphere-Cloud Exchange (PLACE). The PLACE model is a detailed land-surface process model of the partly cloudy atmospheric boundary layer and underlying heterogeneous land surfaces (Wetzel and Boone 1995). Running PLACE requires four index files to specify soil classes, vegetation classes, model constants, and initial soil moisture and temperature. Soil and land cover types must be specified for each of the grid cells in the model domain. Values for the atmospheric variables must be supplied at each time step for each of the grid cells during the length of a simulation. The domain for this study is 50 x 50 grid cells at 1 x 1 km resolution, and it extends both north and south of the Little Washita, OK watershed. Time steps are 300 seconds to match the temporal resolution of atmospheric data. The model period is 9-16 July 1997, encompassing a major precipitation event on 11 July, and pre- and post-event conditions.
In the initial stages of the study, GIS was an excellent tool for organizing and evaluating spatial data required for assimilation into PLACE. GIS themes used in the project include a DEM, land use/land cover grid, and STATSGO polygons. Using Arc/Info's GENERATE function, a polygon coverage of the model domain was developed and unioned to the soil and vegetation coverages. The unioned coverages related each grid cell to a soil type (MUID) and a land cover class. The index files specifying the properties of soil types and land cover in the model domain were developed from the GIS attribute tables for the soils and land cover tables. Soil properties were averages of component and layer properties for each MUID. Vegetation properties were related to the percent cover and classes (i.e. grasses, woody shrubs) of vegetation rather than the species. Extensive processing put the raw atmospheric data from the OK mesonet into a form PLACE could assimilate. A GIS-application was written in the IDL graphics environment to apply the nine widely, irregularly spaced mesonet stations in and around the model domain to a regular grid by Delauney triangulation and then interpolate values for each of the atmospheric variables at each of the model grid cells. Because spline-type interpolation produced negative values, linear interpolation (Cressman weighting) was preferred.
The Future
The single biggest obstacle facing this study is obtaining meaningful results from the PLACE model. As stated earlier, debugging the code to insure it correctly calculates ground heat fluxes is the first priority. Additional problems may arise after fixing the present problem. Subroutines that are called after the ground heat flux subroutine still need to be checked. This obstacle must be overcome before any further improvements are made to the experimental technique.
A number of future improvements could be made to the database. The soil and vegetation properties, particularly the latter, are averages of their MUID or class. County-level and SGP-observed soils data may improve the specification of soil properties because presently, uniform conditions, including elevation changes, are assumed throughout an MUID. Likewise, SGP-observed vegetation data may improve the identification of plant species and quantify their presence in the LW. Further research in the life sciences reference literature may make it possible to derive vegetation properties from averaging species rather than classes in more accurate proportions.
The atmospheric data do not presently take advantage of the USDA micronet. More accurate interpolation of some of the variables would be possible from the denser network within the LW. An interesting comparison would test PLACE results using Doppler radar derived rainfall amounts vs. the interpolated ground observations. Hourly observations of clouds from the National Weather Service would be preferable to estimating the cloud cover fraction used to calculate longwave radiation. In some of the grid cells, forest cover should also be taken into account in the calculation of longwave radiation.
Of course, the effort expended to make improvements to the database should be weighed against any improvement in experiment results. To verify experiment results, both the point gravimetric measurements and remotely sensed soil moisture are required. During SGP, an airborne L-Band (1.4 GHz) microwave sensor took measurements over a narrow swath covering central OK. The images below are of the southwestern portion of the swath.
(Above left) 9 July 1997, pre-storm conditions. The small patch of green and blue is from small shower that took place over the far western portion of the LW. The fringing at the top of the image is from interference from the Oklahoma City airport. (Above right) Early on 11 July 1997 a large storm passed through southern OK. This image is the result of several hours of drying.
Since the point measurements of soil moisture cannot simply be interpolated to a grid, areal maps of soil moisture must be generated using SGP remote sensing data. The point measurements of gravimetric soil moisture and bulk density can be used to relate 1.4 GHz brightness temperature to volumetric soil moisture, so that a proper inversion to volumetric soil moisture can be applied to brightness temperature grids over the model domain. The images above are just that, uncorrected images, not useful for more than display. In September, NASA will release geolocated, noise-corrected grids suitable for quantitative analysis. In the meantime, the images can serve as qualitative checks on model results. GIS can be used to generate and display both model-generated and verification soil moisture grids and create (using ArcView's map calculator) a display of model results. These model results can be quantitatively compared to geography, topography, soil types, and land cover to detect biases in the model-generated grids.
In summary, despite obstacles to conducting these experiments, this study will continue beyond the class project. GIS made possible timely organization of a very large database for assimilation into PLACE. Later, it will make possible spatial analysis of experimental results.
Dingman, S. L., 1996: Physical Hydrology. Prentice Hall, Englewood, New Jersey, pp. 575.
Mahfouf, J.-F, 1991: Analysis of soil moisture from near surface-parameters: A feasibility study. Journal of Applied Meteorology, 30, 1534-1547.
Montieth, J. L., and M. Unsworth, 1990: Principles of Environmental Physics 2nd Ed. Edward Arnold, London, pp. 291.
Rogers, R. R., and M. K. Yau, 1989: A Short Course in Cloud Physics. Pergamon Press, Oxford, UK, pp 293.
Sellers, P. J., W. J. Shuttleworth, J. L. Dorman, A. Dalcher, and J. M. Roberts, 1989: Calibrating the Simple Biosphere model for Amazonian tropical forest using field and remote sensing data. Part I: Average calibration with field data. Journal of Applied Meteorology, 28, 727-759.
Wetzel, P. J., and A. Boone, 1995: A Parameterization for Land-Atmosphere-Cloud Exchange (PLACE): Documentation and testing of a detailed process model of the partly cloudy boundary layer over heterogeneous land. Journal of Climate, 8, 1810-1837.