Advances in Remote Sensing
Vol.2 No.3(2013), Article ID:36009,14 pages DOI:10.4236/ars.2013.23024

On the Use of Landsat-5 TM Satellite for Assimilating Water Temperature Observations in 3D Hydrodynamic Model of Small Inland Reservoir in Midwestern US

Meghna Babbar-Sebens1, Lin Li2, Kaishan Song2, Shuangshuang Xie2

1School of Civil and Construction Engineering, Oregon State University, Corvallis, USA

2Department of Earth Sciences, Indiana University Purdue University Indianapolis, Indianapolis, USA


Copyright © 2013 Meghna Babbar-Sebens et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Received May 11, 2013; revised June 11, 2013; accepted July 11, 2013

Keywords: Data Assimilation; Reservoir Hydrodynamics; Numerical Models; Temperature; Landsat


Accuracy of hydrodynamic and water quality numerical models developed for a specific site is dependent on multiple model parameters and variables whose values are attained via calibration processes and/or expert knowledge. Real time variations in the actual aquatic system at a site necessitate continuous monitoring of the system so that model parameters and variables are regularly updated to reflect accurate conditions. Multiple sources of observations can help adjust the model better by providing benefits of individual monitoring technology within the model updating process. For example, remote sensing data provide a spatially dense dataset of model variables at the surface of a water body, while in-situ monitoring technologies can provide data at multiple depths and at more frequent time intervals than remote sensing technologies. This research aims to present an overview of an integrated modeling and data assimilation framework that combines three-dimensional numerical model with multiple sources of observations to simulate water column temperature in a eutrophic reservoir in central Indiana. A variational data assimilation approach is investigated for incorporating spatially continuous remote sensing temperature observations and spatially discrete in-situ observations to change initial conditions of the numerical model. The results demonstrate the challenges in improving the model performance by incorporating water temperature from multi-spectral remote sensing analysis versus in-situ measurements. For example, at a eutrophic reservoir in Central Indiana where four images of multi-spectral remote sensing data were assimilated in the numerical model, the overall error for the four images reduced from 20.9% (before assimilation) to 15.9% (best alternative after the assimilation). Additionally, best improvements in errors were observed on days closer to the starting time of model’s assimilation time window. However, when the original and updated model results for the water column temperature were compared to the in-situ measurements during the data assimilation period, the error was found to have actually increased from 1.8˚C (before assimilation) to 2.7˚C (after assimilation). Sampling depth differences between remote sensing observations and in-situ measurements, and spatial and temporal sampling of remote sensing observations are considered as possible reasons for this contrary behavior in model performance. The authors recommend that additional research is needed to further examine this behavior.

1. Introduction

In recent years, deterioration of water quality in reservoirs that serve as drinking water sources has become one of the major sources of human health risks. Numerical models have been successfully used to simulate the physical, chemical and biological processes within reservoir systems, and predict the risks of contamination [1-3]. Among all the contributing factors that influence the water quality condition in a reservoir system, water column temperature has significant impacts on the distribution, transportation, and interaction of multiple contaminants such as nutrients, micro-algae, etc. Hence, accurate prediction of water contaminants by numerical water quality models is highly dependent on predictions of water temperature.

Errors in the prediction of temperature can, however, arise from multiple sources, including inaccurate input, inaccurate model parameters, numerical errors during the computation processes, and real-time variations in the system that are not incorporated in the original calibrated model. A data assimilation procedure provides means for integrating real-time observed data from a variety of monitoring sources to improve a model’s prediction accuracy. This improvement is achieved by changing state variables, model inputs, parameter updates, and/or bias correction [4]. In-situ data have been commonly used in the past for assimilation of surface water quality data, including temperature, in numerical models (e.g., [5], etc.). While in-situ monitoring is useful for obtaining multiple temporal observations for multiple depths at a specific X-Y location in the water body, remotely sensed data obtained from satellites provide larger spatial coverage of surface observations than in-situ observations. Hence, a number of data assimilation approaches have been developed to incorporate remote sensing observations for model updates. Data assimilation using remotely-sensed data has been investigated by numerous fields (for example, studies investigating improvements in predictions of soil moisture [6-8], predictions of subsurface soil temperature [9], prediction of snow cover [10] via land surface models, and prediction of hydrodynamic and water quality variables [11-16]. Assimilation of temperature observations has also received considerable interest by numerous researchers (for example, [17-19], etc.). However none of these studies have discussed or investigated how data assimilation using remotely sensed temperature observations affect model accuracy with respect to in-situ temperature observations. This is an important issue for data assimilation procedures since it is very common that information about the errors on the data obtained from a specific data source is unavailable. Additionally in-situ observations may be taken on days and times different from the remote sensed observations, or in-situ observations might not be available at all for the period of interest. This can make it difficult to perform quality control for identifying reliable data, and estimate the error covariance matrix of observation variables. Finally, use of observations with unknown uncertainty may mislead the data assimilation algorithms, as the algorithm tries to “correct” the system state based on these potentially inaccurate observations. It is not known how much of a difference in the model errors would arise when remotesensing data are used to update computationally-expensive numerical models of smaller inland-water bodies.

The specific goals of this study are:

(1) To develop a calibrated 3-Dimensional finite difference numerical model for simulating hydrodynamic processes and temperature in small, inland Eagle Creek Reservoir (ECR) in Central Indiana.

(2) To assimilate water surface temperature retrieved from multi-spectral Landsat-5 TM band 6 images into the hydrodynamic model and adjust the model’s initial conditions via a variational data assimilation approach.

(3) To validate the remote-sensing data-based data assimilation by comparing the updated model results with in-situ observations.

2. Methodology

2.1. Study Area and Data Collection

Eagle Creek Reservoir (ECR) is located about 16 km (10 miles) northwest of Indianapolis, Indiana (Figure 1). It was constructed in 1967 by the city of Indianapolis and was initially used for flood mitigation. A water treatment plant was later constructed and put into service in 1976. The treatment plant takes water directly from the reservoir (approximately 10 million gallons per day) and serves primarily as a source of drinking water supply. ECR is a small and shallow reservoir with normal pool surface area 5.1 km2 and mean depth 5.7 m. It can be separated into three functional areas: the quarry, the northern basin, and the southern basin. The reservoir’s northern and southern basins are separated by a land bridge causeway under 56th street, which allows limited water exchange through an approximately 50 meter opening. The quarry doesn’t have a direct connection with the other two basins in ECR and is considered as an isolate feature. Flow in ECR is supplied by tributaries from the upstream Eagle Creek watershed (426 km2). Four main streams flowing into the reservoir are Eagle Creek, Bush Creek, Fishback Creek and School Branch, with Eagle Creek (mean discharge 4.2m3/s) being the major contributor of the flow.

Seasonal and short-term temperature changes lead to the thermal stratification of the reservoir water. Distinct

Figure 1. Eagle Creek Reservoir morphological conditions (left) and final grid system with 2008 sampling locations (right).

thermoclines that separate the reservoir into warm and cold water zones can prevent the water from mixing. In this case, contaminants tend to accumulate instead of circulating. For example, chlorine can accumulate at the bottom of warm water thermoclines [20]. Also, a lack of rainfall with inadequate mixing of fresh and stagnant water, increased algae growth, deterioration of organic matter as the water warms up, and low wind conditions, have contributed to depletion of DO levels in the reservoir. If a reservoir becomes stratified as a function of temperature, the bottom layer becomes deficient in dissolved oxygen [21].

Bathymetry data for ECR have been measured by Center for Earth and Environmental Sciences (CEES), Indiana University Purdue University Indianapolis (IUPUI). The bottom elevation varies from 223.8 m to 240.5 m above the sea level. Since no existing flow monitoring stations exist on the major tributaries just upstream of the reservoir, inflow discharges into the reservoir from the watershed were obtained from a watershed model—Soil and Water Assessment Tool (SWAT [22])—in this study. The outflow measurements from the dam were obtained from a United States Geological Survey (USGS) gauge Station #03353460 at Clermont (1 km downstream from reservoir). Hourly atmospheric data was obtained from National Climatic Data Center ( Observations from Eagle Creek Airpark/Airport (53842) station were used for this reservoir region. Hourly solar radiation for 2008 was obtained from Indiana State Climate Office ( The closest station-Throckmorton-Purdue Agricultural Center (TPAC)-located in Lafayette, Tippecanoe County, IN, 63 miles Northwest from Indianapolis, was used to obtain solar radiation data. The evaporation in ECR was not known. Hence evaporation rates were based on the measurements of a local water utility company that measured daily evaporation at Carmel, 27 km (17 miles) east of ECR. Data collected by the company indicated an average evaporation value of approximately 5.50 mm/day from June to October, and an average of 4.01 mm/day from November to May. Daily pool elevation data for 2008 was obtained from USGS gauge Station #03353450 in the ECR, just east of the dam.

2.2. Simulation Model

Environmental Fluid Dynamics Code (EFDC) [23,24], a public domain, open source, surface water numerical modeling system for simulating hydrodynamics, and water quality in open-surface water bodies, was used to develop a prediction model for the reservoir. EFDC has been applied to over 100 water bodies to support environmental assessment and management, and regulatory requirements. The EFDC model solves the three-dimensional, vertically hydrostatic, free surface, turbulent averaged equations of fluid with variable density. The model uses a stretched or sigma vertical coordinate and Cartesian or curvilinear, orthogonal horizontal coordinates. The hydrodynamic model also solves dynamically coupled transport equations for turbulent kinetic energy, turbulent length scale, salinity, and temperature.

The 3-D stretched sigma grids implemented with the EFDC use the following transforming function to calculate an adjusted vertical coordinate, Z, from the bottom elevation and water surface elevation:


In Equation (1), Z* = original physical vertical coordinate, h = bottom elevation and ζ = water surface elevation. After the physical vertical coordinate system is stretched, the total depth is evenly distributed into equal depths of individual layers, for all the X-Y grid locations within the research domain.

The continuity equation used in EFDC is given by Equation (2), in which H = water depth, u and v = horizontal velocity components in x and y direction respectively, w = vertical velocity component in z direction; QH = the volumetric source and sink term concerning rainfall, evaporation and infiltration. The conservation of momentum equations are given in Equations (3) and (4), in which f = Coriolis factor, p = the water column hydrostatic pressure; patm = the kinematic atmospheric pressure; Av = vertical turbulent momentum diffusion coefficients, and Qu and Qv = momentum source-sink terms.




In the transport equations for salinity and temperature (Equations (5) and (6)) the source and sink terms are given by QS and QT, which consist of sub-grid scale horizontal diffusion and thermal sources and sinks, and Aw is the vertical turbulent diffusivity. 



The water surface and bed boundary conditions for heat transport are given by Equations 7 and 8. For water surface:


For bed:


Short wave solar radiation at the bed is defined as:


where Jb = net long-wave back radiation; Jc = convective heat transfer; Je = evaporation heat transfer; Cpw = specific heat of water; Hb = active thermal thickness of the bed; Tb = bed temperature; Ib = short-wave solar radiation at the bed; ρb = bed density; Cpb = specific heat of the water-solid bed mixture; Chb = dimensionless convective heat exchange coefficient; Tbl = bottom layer water temperature; Is = solar radiation at the water surface, r = distribution factor; βf and βs = fast and slow-scale attenuation coefficients (Caliskan, 2008).

Eight Equations (2)-(9) provide a closed system for the variables u, v, w, p, z, r, S, and T. The vertical turbulent viscosity and diffusivity and the source and sink terms are also specified [23]. EFDC uses mode-splitting [25] for separating the vertically integrated equations (external free surface gravity wave or barotropic mode) from the vertical structure equations (internal shear or baroclinic mode). It calculates free surface elevation by solving the velocity transport separately from the 3D calculation of velocity and thermodynamic properties.

2.2.1. EFDC Grid Generation

Multiple grid sizes and time steps for representing the physical system were explored in the study to assure the accuracy of model results as well as the efficiency of the model. The Courant-Friedrichs-levy condition (CFL condition), a necessary condition for convergence while solving certain partial differential equations (usually hyperbolic PDEs) using a two-time-level numerical scheme, was used to test the suitability of grid sizes and time steps. CFL condition is expressed as Equation (10):


where U = velocity, ∆t = time step, ∆x = cell size. The necessary restriction for grid size and time step to ensure numerical convergence and stability is γ < 1.

The final grid setup (Figure 1) for the numerical model consisted of expanding grids with minimum grid size 40 m to maximum 60 m. The expanding factor of 1.005 was chosen to expand grid sizes from the focal point, which was the water intake. This location was chosen to accurately simulate the most complex flow condition happening close to the causeway under the 56th Street land bridge, and the drinking water intake. A total of 2401 grid cells were developed to represent the physical domain in the modeling domain. This grid best represented the shape of the reservoir shoreline compared to other grids. Since EFDC only recognizes flow through cell faces, cells connecting with each other by corners do not exchange any mass or momentum across the corners. For this grid system, the time step of the finite difference model was set up to be two seconds considering the model stability and the computational burden.

2.2.2. Initial/Boundary Conditions

The depth of water through the reservoir was measured by CEES (Center for Earth and Environmental Science, IUPUI) with sonar equipments. Bottom elevation was based on the bathymetry data, and the initial water surface elevation (or pool elevation (PE)) for the model was chosen to be 240.56 m according to USGS (station near the dam) measurement on January 1st, 2008. A uniform initial water temperature of 3.4˚C throughout the reservoir was assumed for the model, based on the measurements at Mill Creek USGS gauge station near Manhattan, IN on the same day.

Hydrodynamic boundary conditions of ECR included 1) time-series inflow discharge from twelve tributaries of Eagle Creek simulated by the Soil and Water Assessment Tool - based watershed model (Figure 2), 2) outflow discharge through the water intake and ECR dam, 3) wind speed and direction, and 4) atmospheric data including precipitation and evaporation. The climatic data for 2008 was obtained from National Climate Data Center (NCDC) and Indiana State Climate Office. Eagle

Figure 2. Eagle creek watershed tributaries that interface with the numerical modeling grid of ECR.

Creek Airpark station (station ID 53842) was used to obtain data on wet/dry air bulb temperature, precipitation, relative humidity, air pressure, wind speed and direction, and cloud cover. All data were appropriately formatted and imported into the EFDC modeling system. Figure 2 also shows the boundary where the watershed tributaries join the EFDC modeling grid for ECR.

2.3. Remote Sensing Observations

One of the primary objectives of this study was to examine the advantages of using remote sensing observations for updating the hydrodynamic numerical model of inland freshwater water body, hence, multispectral data were obtained from the Landsat 5 satellite. Landsat 5 Thematic Mapper (TM) has been previously used to derive water temperature with measurable success. For example, Lathrop and Lillesand [26] analyzed the utility of Landsat-5 TM thermal IR (band 6) for measuring and mapping water temperature of the Great Lakes, and concluded that the TM derived water surface temperature has a root mean square error of less than 1˚C. Wouthuyzen et al. [27] used 13 dates of Landsat-5 TM thermal images to investigate the seasonal variation of temperature in the Omura Bay in western Kyushu Island of Japan, and the results from this study indicates that water surface temperature could be estimated at an accuracy of 0.551˚C, 0.371˚C, 0.351˚C, and 0.331˚C for spring, summer, autumn, and winter, respectively. Schneider and Mauser [28] found that lake surface temperature could be estimated at an accuracy of 0.55˚K. Schneider et al. [29] have used Landsat TM thermal band to estimate water surface temperature of Lake Constance and achieved an accuracy of 0.53˚K.

In this research, the spatial resolution for TM band 6 images was re-sampled from 120 m to 30 m spatial resolution. Four of the TM band 6 thermal images obtained from the satellite (on dates Aug 7, Aug 23, Sept. 24, and Oct. 10 2008, and at local times 16:09, 16:08, 16:07, and 16:07 hours, respectively) were used to convert spectral radiance to water surface temperature based on the Planck’s law (Equation (11)). In Planck’s law, the long wave radiation emitted from the land surface was in proportion to its temperature as:


where 𝜆 is the wave length; K1 and K2 are the calibration constants as 607.76 Watts/(m2·ster·μm) and 1260.56 Kelvin respectively [21,30]; L is the spectral radiance in watts/(m2·ster·μm). These pre-launch calibration constants from the empirical models are used under the assumption that the down-welling radiance and atmospheric transmissivity are constant in space throughout the study area, and thus are applied to calculate surface temperature for each image pixel [31]. Studies using this method found the root mean square error (RMSE) in water surface temperature to be less than 1K of retrieving land surface temperature [32]. This error is even smaller for surface water systems because of their better homogeneity in surface temperature.

2.4. Data Assimilation Algorithms

Existing methods of data assimilation are based on two types of approaches for finding the best estimates of state variables, input variables and boundary conditions from (noisy) observations given a (noisy) model [33]. The first approach uses a “direct observer” and provides a fourdimensional data assimilation scheme, whereas the second approach uses a “dynamic observer” and provides a sequential data assimilation scheme. The commonly used direct observer data assimilation approaches include Direct Insertion [34,35], Statistical Correction [34], Nudging [34,36], and Kalman Filter related methods. These “direct observer” methods adjust the model by continuous update according to the observation in the previous time step. The “dynamic observer” data assimilation approaches adjust the state variables at the beginning of each assimilation window so that model predictions over that time period correspond with the observations. “Dynamic observer” techniques are especially useful for problems that are driven by accuracy of initial conditions. “Dynamic observer” techniques can be posed as optimization problems with strong constraints (variational methods) or weak constraints (dual variational or representer methods). Dente et al. [37] successfully used a variational method to assimilate ASAR and MERIS satellite data into a wheat model and improved the wheat yield mapping. They employed a cost function to measure the error between model outputs and observations. In order to search for the optimal configuration of the model’s initial conditions, they minimized the cost function in an optimization method by varying the initial conditions. The optimization method constrained the variability of initial conditions to lie between expected values of initial conditions. The model was reinitialized within this range of initial conditions and the optimum model simulation of wheat yield was obtained. In the research done by Ines et al. [38], remote sensing data from two Landsat-7 Enhanced Thermatic Mapper Plus (ETM+) band 6 images were assimilated into a soil-water-atmosphere-plant model using variational data assimilation method. A genetic algorithm (GA) was used in data assimilation to modify model initial conditions and for optimization of water management strategies. A variational data assimilation method has also been used to assimilate a sequence of satellite images into a simple transport-diffusion model to simulate the ocean surface current [39]. Although data assimilation algorithms have been widely applied in environmental studies, few studies related to surface water systems [14,39] have used “dynamic observer” data assimilation approaches into the hydrodynamic models.

The main objective of data assimilation approach in this study was to change the values of initial conditions of temperature in the reservoir at the beginning of the simulation, so that the cumulative error between model predictions of temperature and observations at a future time was minimized. A cost function that estimates the error between predictions and observations in the surface layer, at specific time periods and at specific locations, was constructed and used to assist the optimization algorithm (based on a genetic algorithm) in modifying the assumed initial conditions of temperature in the various layers of the reservoir. The optimization algorithm was run for multiple iterations to improve the initial conditions.

To assimilate data obtained from the Landsat-5 TM+ images, 300 random locations were identified in the reservoir (Figure 3), where the model results and optical observations were used to evaluate the model performance. The finite difference model computed and generated results in 2401 cells distributed over five vertical layers within ECR. Since the remotely sensed observations are most applicable to the water surface conditions, the very top layer of the finite difference grids was chosen for comparing model predictions to the corresponding satellite-derived temperature observations.

A variational data assimilation algorithm (Figure 4) was used in this study to incorporate the remote sensing observations within the simulation model. The data assimilation window was defined as the time period from the initial condition to the last observation time (Figure 4). The cost function was estimated using the relative root mean square error (RRMSE) calculated from both remote sensing outputs and model outputs in the topmost layer of all 300 random locations (Equation (12)). RRMSE calculated for each of the four observation days was then

Figure 3. Random locations identified in the reservoir with TM image on Aug 7, 2008.

Figure 4. Variational data assimilation of temperature observations.

equally weighted and summed to obtain an overall RRMSE for the cost function (Equation (13)).



where Xo,i,j is value of the field observed model parameter at the ith location and jth day; Xm,i,j is the model output at the ith location and jth day; n is the number of observations at the the various locations (i.e. 300 in this case study); is the mean value of all the spatial observations; j is the number of days when remote sensing observations were acquired.

An evolutionary optimization scheme based on the single objective genetic algorithm (GA [40]) was used in this study for minimizing the cost function (Equation (12)). The five decision variables were defined as the percentage change in the initial conditions of temperature for each of the five vertical layers. Nine uniformly distributed values of percent change were chosen between the range of minimum and maximum values of percent change. The first remote sensing image was available on the 219th day (August 7th) of 2008. Hence, the data assimilation window was started from 213th (August 1st) day to allow a warm-up period for EFDC. Since the hydrodynamic numerical model is computationally timeconsuming, use of evolutionary optimization-based data assimilation methods can pose additional computational challenges in effectively searching for alternatives. Hence, we limited the GA population size to a small value of 16, and set the maximum number of generations to 10. However, small population size and generation size can further limit the search efficiency. Hence, we conducted the optimization scheme over multiple stepby-step experiments. In the first GA experiment, the initial range of percentage change in water temperature for all the five vertical layers was allowed to vary from −20% (minimum percent change) to 20% (maximum percent change) of the default initial temperature on 213th day. In the following experiments, this range was individually modified for each vertical layer to allow search for better values of percent change for each of the layers near the values of the best alternative found in the previous experiment. In total five experiments were implemented to the GA algorithm (Table 1). The crossover and mutation rate were set to 0.90 and 0.05, respectively.

3. Results and Discussion

3.1. Hydrodynamic Model Calibration

3.1.1. Water Surface Elevation and Mass Balance

The period January 1st, 2008 to July 31st, 2008 were used for the calibration process, just before the first set of remote sensing observations was available. The initial root mean square error (RMSE) for the calculated pool elevation was estimated to be 1.036 meters. These discrepancies arose due to errors in the estimated tributaries inflows (obtained from SWAT hydrologic model), which indicated need for further adjustment of inflows in order to better match the observed pool elevation. Also, since such large discrepancies in water surface elevation could lead to significant mismatches in water column temperature vertical profiles, a further adjustment in SWAT model outputs of flow discharges was made. The first step in adjusting the SWAT inflows used for the EFDC model included calculation of daily net storage in the reservoir from the inflow and outflows based on the following mass balance Equation (14):


Table 1. Decision variables ranges (in %) for the five vertical layers in the five GA experiments (“Exp.”). Choice of ranges for 2nd, 3rd, and later experiments is dependent and based on the initial condition of best alternative found in the previous experiment.

Once the daily net storage was estimated, it was added to the volume of the water in the reservoir to estimate the volume of the reservoir at the beginning of the next day (Equation (15)):


A rating curve was then developed using the bathymetry data, to estimate relationships between reservoir volume and pool elevation. This rating curve provided means to estimate the reservoir volume at any day based on the pool elevation measured on that day. The difference between the calculated and measured water volume in the reservoir was weighted according to the tributary discharge for that day, and was subtracted from the daily tributary inflows (Equation (16)).

After tributary inflows were adjusted, the new net stor-


age of water within the reservoir and the corrected total volume on day i were calculated. At the beginning of day i + 1, the initial volume was the summation of corrected water volume in the day i and the net storage from inflow and outflow. The total difference between measured and calculated water volume was again used in the correction step, generating the corresponding tributary flows for day i + 1. This adjustment was applied to all daily inflows until the end of the model simulation period. The schematic procedure of computation is shown in Figure 5. Though this adjustment produced the water surface pool elevation identical to the measured data, it also resulted in some negative values in the adjusted tributary flows. In order to correct the negative tributary flows, all the negative flows were replaced by a very small positive discharge of 0.00001m3/s. This resulted in a small discrepancy (RMSE of 0.07m) in pool elevation estimated from the rating curve. This error after flow correction was within the reasonable range in scale.

The adjusted and corrected SWAT tributary inflows were then used as inputs into the reservoir model (EFDC). Figure 6 compares the EFDC simulated water surface elevation with the measured pool elevation. Since EFDC simulates at much smaller time step than the daily timestep of the SWAT model, the EFDC output at the end of each day was taken as the corresponding simulated value for the water surface elevation on that day in Figure 6. Comparison of the modeled results and the field measurements during the first seven months in 2008 produced a root mean square error (RMSE) of 0.029 meter, which was 0.171 meter below the 0.200 meter instrumental accuracy of the Acoustic Doppler Current Profiler used to measure the bathymetry (Lobligeois, 2009).

3.1.2. Water Column Temperature

Other than water surface elevation, water column temperature was calibrated in full hydrodynamics. The insitu measurements of 2008 were collected at 54 different (X, Y) locations throughout the reservoir, and with multiple measurements at several depths in each location. The YSI probes that were used to measure the water column temperature at these sampling locations had an

Figure 5. Water surface elevation as a result of flow adjustment (The bottom figure is an enlarged sub-section of the top figure).

instrumental error of ±0.15˚C (YSI, 2011). These measurements were used for model calibration and validation. During the model calibration period, a series of calibration parameters were adjusted using EFDC to achieve the best agreement with the in-situ measurements of water column temperature (Table 2). In this research, water temperature calibration used 23 sampling events between May 22, 2008 and July 30, 2008.

The temperature calibration results for 16 measurement events are shown in Figure 7. The densely dashed horizontal lines show the water surface elevation of each corresponding water column; the solid greenish-brown colored horizontal lines represent the bottom elevation. It indicates very good correlation between the measured (blue curves with square markers) and the modeled (brown smoothed curves with no markers) for water column temperature. Statistical analysis shows an overall RMSE of 1.279˚C, which is within the satisfaction range for water temperature simulation using EFDC. From the statistical view, the water temperature simulation has the most significant percentage error at monitoring Station ECRAT_A4. This phenomenon can be explained because this station is located at the upstream riverine area of ECR, and at the downstream of the confluence of an unnamed tributary with the reservoir. In summer, it is assumed that water around this shallow region of the reservoir has higher temperature than nearby tributary flows from the watershed. The tributary flows bring cooler wa-

Table 2. Water temperature calibration parameters.

Figure 6. Comparison between measured pool elevation (PE), EFDC simulated pool elevation, and pool elevation calculated from daily mass balance based on SWAT daily flows and rating curves.

Figure 7. Calibrated temperature vertical profiles for ECR.

ter from the watershed into this area, where more stable and warmer water conditions have been created. However, due to the lack of data from the watershed model, the reservoir boundary condition for water temperature was determined based on the USGS gauge station #03354000 at White River near Centerton, IN (approximately 25 miles south from ECR). Therefore, the modeled temperature has high deviation from measured temperature at this station. With the calibrated water column temperature provided by EFDC calibrated water column temperature, a remote sensing data assimilation approach was then implemented during the simulation time period from August 1st to October 10th, 2008.

3.2. Remote Sensing Data Assimilation

Since one of the main objectives of this study was to investigate how a model updated via assimilation of remotely sensed observations of temperature would perform with respect to in-situ observations of temperature, the optimization algorithm calculated the error cost function (Equation (13)) based on four remote sensing images during 2008 (August 7th, August 23rd, September 24th, October 10th). At this point, it is important to note that the in-situ observations were not available on the days the remote sensing observations were available, which is a practical monitoring scenario that can arise. Hence errors were calculated for different set of days for the two types of observations.

During the remote sensing data assimilation processes, EFDC was set up in a restart mode with the restarting time at August 1st, which was the 213th day of the simulation time period. The EFDC model was reinitialized with the water temperature initial condition adjusted at 213rd day of 2008. During this process of reducing the variational range of the decision variables, it was discovered that the decreasing water temperature values in initial condition provide better results, which were closer to the remote sensing observations. The results in assimilating water temperature data from the TM images successfully reduced the overall error fitness function (Equation (13)) from 20.9% (before assimilation) to 15.9% (best alternative after the assimilation). The best alternative had decision variable values that indicated −55%, −41.25%, −60%, −60%, and −51.25% reductions in the initial temperatures of layers 1 (deepest), 2, 3, 4, and 5 (surface layer). Also, the improvements in model performance reduced rapidly with the progression of time, within the assimilation window. For example, for 7th Aug 2008, the best model’s relative RMSE improved from 22.0% to 12.1% (i.e. 9.9% relative RME reduction), whereas for 10th Oct 2008 the best model’s relative RMSE improved from 17.3% to only 14.8% (i.e., 2.5% relative RMSE reduction).

Figure 8 shows comparison between the best modeled surface layer water temperature at 300 random locations in the reservoir (circles) and the value retrieved from Landsat TM images, on the four observation days in 2008. The temperature outputs after data assimilation application showed better consistency with remote sensed observations. These results also indicate a reasonably fair performance of the genetic algorithm optimization scheme in the variational data assimilation process, in spite of the use of smaller population sizes and fewer generations because of computational burden of the numerical model.

The errors were next calculated between modeled and in-situ measurements for the model before data assimilation and for the best model after data assimilation, and for the time period covering the data assimilation window. Table 3 lists these comparisons of at various in-situ stations for the entire water column. For 15 out of 23 X-Y locations (i.e., 65% of the station locations), the water column error unexpectedly increased after the data assimilation. The overall average RMSE of calibrated EFDC model (i.e., the model before data assimilation) in comparison with the in-situ observations, for the water column and during the data assimilation window, was estimated to be 1.8˚C. After data assimilation, this overall average RMSE of the best updated EFDC model increased to 2.7˚C. This indicated that even though the data assimilation improved the model predictions of surface layer temperatures with respect to the remote sensing temperature observations, the updated model’s performance with respect to the in-situ observations in the water column actually worsened by 0.9˚C. This conflicting change in errors indicates that data assimilation based on the remote sensing derived water surface temperature can

Figure 8. Water temperature data assimilation results: Circles-EFDC generated temp at 300 points without DA adjust; triangles-DA adjusted tempt at 300 points.

Table 3. Statistical summary of remote sensing data assimilation (RSDA) on water column temperature.

force the model farther away from the in-situ measurements of water column temperatures.

3.3. Comparison between Remote Sensing Data and in-Situ Measurements

The preliminary data assimilation results proved that remote sensing data can force model updates in directions that have conflicting model performance with respect to in-situ data. These differences might come from random or unknown human error, which are not easily estimated, or from sampling biases. Hence, it was decided to investigate the characteristics of available water surface observations to measure any potential bias in the measurements. Since, the two sets of data were taken on different days and none of the in-situ observations accurately overlapped the spatial locations of previous in-situ observations (i.e., one (x,y) location was sampled only once in 2008), it was decided to use a grouping scheme to allow comparisons between measurements. In-situ observations were clustered into four representative physical regions (Regions A, B, C, and D) identified within the reservoir, which provided means to create a time series of in-situ observations at a specific region. Each in-situ observation taken on a specific day and at the shallowest depth (0.25 meters) of a specific (x, y) location was used to represent the in-situ-based surface temperature observation of the relevant region. Similarly, the four temporal remotely-sensed observations in the respective regions were also used to estimate the relevant remote sensing derived water surface temperature (i.e., skin temperature representative of depth less than 0.05 meter) for the various regions over the four days monitored. This was done by taking the average of all remote sensing derived temperature observations at the (x, y) locations of the in-situ observations in the relevant region. Region A and B were located in the northern basin of the reservoir, and Region C and D were located in the southern basin (Figure 9). Figures 10 (a)-(d) compare the available observations of four remote sensing derived

Figure 9. Four geographical sampling regions in ECR.

Figure 10. (a)-(d) Comparison between CEES measurements and remotely derived temperature.

temperature with the in-situ temperature observations for the four regions. The instrument error of ±1˚C interval for remote sensing data and ±0.15˚C interval for in-situ measurements are also included as lower and upper bars in these figures. It can be seen that in-situ observations are generally higher than remote sensing derived observations for all the four regions, thereby, showing a bias in the remote-sensing observations which are generally cooler than in-situ measurements. Statistical results indicated an average RMSE of 3.33˚C between in-situ observations and remote-sensing observations, with in-situ observations having higher temperatures. Since synchronous measurements from the satellite and in-situ sensors are not available, we compared the water temperature measurements with the dry bulb temperature of air. Figure 10(a) compares these water temperatures in region A with the air temperature. From this comparison it can be assessed that all four remote sensing observations were taken on days with air temperature mostly cooler than air temperatures during the neighboring days when in-situ observations were taken. The same can be stated about regions B, C, and D, indicating that since all four days when remote sensing data was available, the air temperature was also cooler. The lack of enough remote sensing data during the warmer days could have biased the data assimilation towards cooler temperatures, and hence increasing the errors with respect to warmer in-situ temperatures.

4. Conclusions

This research explored the limitations in using remote sensing data for data assimilation in a finite difference hydrodynamic model of a small, inland freshwater reservoir. A genetic algorithm-based variational data assimilation approach successfully generated model initial conditions which reduced the difference between the model outputs and the remote sensing observations. However, the adjusted/updated model produced predicttions of temperature that had higher errors with respect to the in-situ measurements. Therefore, the assumption that the data assimilation with remote sensing temperature observations will also always improve model performance with respect to field measurements was rejected. Following are the challenges that could potentially be the cause behind the observed discrepancies:

1) The data assimilation algorithm was performed based on model outputs and remote sensed TM images. Horizontally, it was assumed that each random location was fully representative of the area round within one 30 m remote sensing pixel and model grid. Since the spatial resolution is smaller for TM images than the model grid system, the unit pixel for each data set would be mismatched in spatial scale. For example, when two data sets were overlapped with each other, there would be more than one pixel lain over the single model grid. When the random locations were picked up by points, the corresponding value extracted from the TM images would not necessarily be the dominant pixel value located within its model grid size. Similarly, when identifying the in-situ measurement locations, the corresponding values would not necessarily be dominant pixel value located within its model grid size from the TM image. Thus assimilating corresponding remote sensed water temperature into the model grids would introduce errors during the horizontal integrating process. 

2) When comparing the remote sensing data with insitu measurements, it was assumed that the measurements at 0.25 meter depth would be comparable to remote sensed water surface temperature. However, as the remote sensing technology can only represent the top surface of the water column, the difference between the water temperature on the skin and that at 0.25 meter depth was assumed to be insignificant. However, based on our results, this difference warrants additional research and should be examined in relation to the temperature gradient at the time of in-situ measurements to avoid the influence cause by water depth.

3) We also observed that the four remote sensing observations were generally taken on cooler days, in comparison to the in-situ observations on neighboring days. The effect of lack of additional remote sensing data during the warmer days needs to be also investigated in future research.

We propose that since currently there is no method for accurately measuring the systematic and random human errors in either remote sensing or in-situ measurements, both remote sensing and in-situ measurements need to be assimilated at the same time during the data assimilation period. Future studies are also needed to investigate (a) how a multi-objective type approach can be applied for variational data assimilation of data from multiple sources, and (b) how discrepancies in spatial and temporal scales of remote-sensing observation, in-situ observations, and model discretization needs to be resolved when they are used in combination in data assimilation approaches.

5. Acknowledgements

The authors would like to thank the funding agencies, Central Indiana Water Resources Partnership and National Aeronautics and Space Administration (Award ID# NNX10AH13G), for supporting this research.


  1. United States Environmental Protection Agency (USEPA), “Total Maximum Daily Load for Dissolved Oxygen and Nutrients to Mashapaug Pond,” Rhode Island, 2002.
  2. K. Jin and Z. Ji, “Case Study: Modeling of Sediment Transport and Wind-Wave Impact in Lake Okeechobee,” Journal of Hydraulic Engineering, Vol. 130, No. 11, 2004, pp. 1055-1067. doi:10.1061/(ASCE)0733-9429(2004)130:11(1055)
  3. Tetra Tech, “Hydrodynamic and Water Quality Modeling Report for Lake Lanier,” Georgia, 2009.
  4. H. Moradkhani, “Hydrologic Remote Sensing and Land Surface Data Assimilation,” Sensors, Vol. 8, No. 5, 2008, pp. 2986-3004. doi:10.3390/s8052986
  5. J.Q. Mao, J. H. Lee and K. W. Choi, “The Extended Kalman Filter for Forecast of Algal Bloom Dynamics,” Water Research, Vol. 43, No. 17, 2009, pp. 4214-4224. doi:10.1016/j.watres.2009.06.012
  6. X. Li, T. Koike and M. Pathmathevan, “A Very Fast Simulated Re-Annealing (VFSA) Approach for Land Data Assimilation,” Computers & Geosciences, Vol. 30, No. 3, 2004, pp. 239-248. doi:10.1016/j.cageo.2003.11.002
  7. S. Boussetta, T. Kolke, K. Yang, T. Graf and M. Pathmathevan, “Development of a Coupled Land-Atmosphere Satellite Data Assimilation System for Improved Local Atmospheric Simulations,” Remote Sensing of Environment, Vol. 112, No. 3, 2008, pp. 720-734. doi:10.1016/j.rse.2007.06.002
  8. L. Zhu, J. M. Chen, Q. Qin, J. Li and L. Wang, “Optimization of Ecosystem Model Parameters Using SpatioTemporal Soil Moisture Information,” Ecological Modelling, Vol. 220, No. 18, 2009, pp. 2121-2136. doi:10.1016/j.ecolmodel.2009.04.042
  9. A. Olioso, H. Chauki, D. Courault and J. Wigneron, “Estimation of Evapotranspiration and Photosynthesis by Assimilation of Remote Sensing Data into SVAT Models,” Remote Sensing of Environment, Vol. 68, No. 3, 1999, pp. 341-356. doi:10.1016/S0034-4257(98)00121-7
  10. A. G. Slater and M. P. Clark, “Snow Data Assimilation via an Ensemble Kalman Filter,” Journal of Hydrometeorology Vol. 7, No. 3, 2006, pp. 478-493. doi:10.1175/JHM505.1
  11. J. Yang and F. X. LeDimet, “Variational Data Assimilation in the Transport of Sediment in River,” Science in China Series D-Earth Sciences, Vol. 41, No. 5, 1998, pp. 473-485. doi:10.1007/BF02877738
  12. D. Seo, V. Koren and N. Cajina, “Real-Time Variational Assimilation of Hydrologic and Hydrometeorological Data into Operational Hydrologic Forecasting,” Journal of Hydrometeorology, Vol. 4, No. 3, 2003, pp. 627-641. doi:10.1175/1525-7541(2003)004<0627:RVAOHA>2.0.CO;2
  13. M. N. Madsen, J. N. Hartnack, C. Skotner and A. Berg, “Water Quality Surveillance and Early Warning in Surface Waters—Intergration of Mathematical Models and On-line Monitoring,” Proceedings of the 7th International Conference on Hydroinformatics, Nice, 4-8 September 2006, pp. 1-23.
  14. G. Panteleev, A. Proshutinsky, M. Kulakov, D. A. Nechaev and W. Maslowski, “Investigation of the Summer Kara Sea Circulation Employing a Variational Data Assimilation Technique,” Journal of Geophysical ResearchOceans, Vol. 112, No. C4, 2007, Article ID: C04S08.
  15. A. Voutilainen, T. Pyhalahti, K.Y. Kallio, J. Pulliainen, H. Haario and J. P. Kaipio, “A Filtering Approach for Estimating Lake Water Quality from Remote Sensing Data,” International Journal of Applied Earth Observation and Geoinformation Vol. 9, No. 1, 2007, pp. 50-64. doi:10.1016/j.jag.2006.07.001
  16. Alabama Department of Environmental Management, “Water Quality Branch, Final Coosa River Basin Total Maximum Daily Loads for Nutrients, OE/DO and PH,” Final Report, 2008.
  17. T. Ezer and G. L. Mellor, “Simulations of the Atlantic Ocean with a Free Surface Sigma Coordinate Ocean Model,” Journal of Geophysical Research, Vol. 102, 1997, pp. 15647-15657. doi:10.1029/97JC00984
  18. C. L. Keppenne and M. M. Rienecker, “Assimilation of Temperature into an Isopycnal Ocean General Circulation Model Using a Parallel Ensemble Kalman Filter,” Journal of Marine Systems, 2003, Vol. 40-41, pp. 363-380. doi:10.1016/S0924-7963(03)00025-3
  19. A. Troccoli and K. Haines, “Use of the TemperatureSalinity Relation in a Data Assimilation Context,” Journal of Atmospheric and Oceanic Technology, Vol. 16, No. 12, 1999, pp. 2011-2025. doi:10.1175/1520-0426(1999)016<2011:UOTTSR>2.0.CO;2
  20. United States Environmental Protection Agency (USEPA), “Uncovered Finshed Water Reservoirs,” EPA Guidance Manual, 1999.
  21. D. Goodin, “Mapping the Surface Radiation Budget and Net Radiation in a Sand Hills Wetland Using a Combined Modeling/Remote Sensing Method and Landsat Thematic Mapper Imagery,” Geocarto International, Vol. 10, No. 2, 1995, pp. 19-29. doi:10.1080/10106049509354488
  22. J. G. Arnold, R. Srinivasan, R. S. Muttiah and J. R. Williams, “Large Hydrologic Modeling and Assessment: Part I. Model Development,” Journal of American Water Resources Association, Vol. 34, No. 1, 1998, pp. 73-89. doi:10.1111/j.1752-1688.1998.tb05961.x
  23. J. M. Hamrick, “A Three-Dimensional Environmental Fluid Dynamics Computer Code: Theoretical And Computational Aspects,” Applied Marine Science and Ocean Engineering, Virginia Institute of Marine Science School of Marine Science, The College of William and Mary, Gloucester Point, Special Report No. 317, 1992.
  24. United States Environmental Protection Agency (USEPA), Environmental Fluid Dynamics Code (EFDC), 2010.
  25. T. J. Simons, “Verification of Numerical Models of Lake Ontario. Part I. Circulation in Spring and Early Summer,” Journal of Physical Oceanography, Vol. 4, 1974, pp. 507- 523. doi:10.1175/1520-0485(1974)004<0507:VONMOL>2.0.CO;2
  26. R. G. Lathrop Jr. and T. M. Lillesand, “Calibration of Thematic Mapper Thermal Data for Water Surface Temperature Mapping: Case Study on the Great Lakes,” Remote Sensing of Environment, Vol. 22, No. 2, 1987, pp. 297-307. doi:10.1016/0034-4257(87)90063-0
  27. S. Wouthuyzen, K. Gotoh and S. Uno, “Seasonal Sea Surface Temperature of the Omura Bay, Japan Estimated from Multi-Date Landsat-5 TM Thermal Infrared Band,” Proceedings of 1993 International Geoscience and Remote Sensing Symposium, Tokyo, 18-21 August 1993, pp. 146-148.
  28. K. Schneider and W. Mauser, “Processing and Accuracy of Landsat Thematic Mapper Data for Lake Surface Temperature Measurement,” International Journal of Remote Sensing, Vol. 17, No. 1, 1996, pp. 2029-2041. doi:10.1080/01431169608948757
  29. K. Schneider, H. Bach and W. Mauser, “Determination of Water Quality Parameters in Lake Constance Using LANDSAT-TM Images,” SPIE Proceedings of Remote Sensing of Vegetation and Sea, Taormina, 17 January 1997, pp. 149-158.
  30. National Aeronautics and Space Administration (NASA), Landsat 7 Science Data Users Handbook, 2009, pp. 106- 128.
  31. C. A. Laymon and D. A. Quattrochi, “Estimating Spatially Distributed Surface Fluxes in a Semi-Arid Great Basin Desert Using Landsat TM Thermal Data,”In: D. A. Quattrochi and J. C. Luvall, Eds., Thermal Remote Sensing in Land Surface Processing, CRC Press, Boca Raton, 2004, pp. 133-159.
  32. J. A. Sobrino, J. C. Jimenez-Munoz and L. Paolini, “Land Surface Temperature Retrieval from LANDSAT TM 5,” Remote Sensing of Environment, Vol. 90, No. 4, 2004, pp. 434-440. doi:10.1016/j.rse.2004.02.003
  33. J. P. Walker and P. R. Haouser, “Hydrologic Data Assimilation,” In: A. Aswathanarayana, Ed., Advances in Water Science Methodologies, A.A. Balkema, The Netherlands, 2005, 230 p.
  34. P. R. Houser, W. J. Shuttleworth, J. S. Famiglietti, H. V. Gupta, K. H. Syed and D. C. Goodrich, “Integration of Soil Moisture Remote Sensing and Hydrologic Modeling Using Data Assimilation,” Water Resources Research, Vol. 34, No. 12, 1998, pp. 3405-3420. doi:10.1029/1998WR900001
  35. A. R. Robinson and P. F. J. Lermusiaux, “Overview of Data Assimilation,” Harvard Reports in Physical/Interdisciplinary Ocean Science, 2000.
  36. P. Malanotte-Rizzoli and R. E. Young, “Assimilation of Global versus Local Data Sets into a Regional Model of the Gulf Stream System .1. Data Effectiveness,” Journal of Geophysical Research-Oceans Vol. 100, No. C12, 1995, pp. 24773-24796. doi:10.1029/95JC02580
  37. L. Dente, G. Satalino, F. Mattia and M. Rinaldi, “Assimilation of Leaf Area Index Derived from ASAR and MERIS Data into CERES-Wheat Model to Map Wheat Yield,” Remote Sensing of Environment, Vol. 112, No. 4, 2008, pp. 1395-1407. doi:10.1016/j.rse.2007.05.023
  38. A. V. M. Ines, K. Honda, A. D. Gupta, P. Droogers and R. S. Clemente, “Combining Remote Sensing-Simulation Modeling and Genetic Algorithm Optimization to Explore Water Management Options in Irrigated Agriculture,” Agricultural Water Management, Vol. 83, No. 3, 2006, pp. 221-232. doi:10.1016/j.agwat.2005.12.006
  39. G. K. Korotaev, E. Huot, F.-X. Le Dimet, I. Herlin, S. V. Stanichny, D. M. Solovyev and L. Wu, “Retrieving Ocean Surface Current by 4-D Variational Assimilation of Sea Surface Temperature Images,” Remote Sensing of Environment, Vol. 112, No. 4, 2008, pp. 464-1475. doi:10.1016/j.rse.2007.04.020
  40. D. E. Goldberg, “Genetic Algorithms in Search and Optimization and Machine Learning,” Addison-Wesley, Boston, 1989.