The Modified SEBAL for Mapping Daily Spatial Evapotranspiration of South Korea Using Three Flux Towers and Terra MODIS Data

Evapotranspiration (ET) is expected to increase by a considerable amount because of the impact of future temperature increase. Nowadays, the daily to seasonal ET maps can be used to provide information for a sustainable and adaptive watershed eco-environment. This study attempts to estimate the spatial ET of South Korea (99,900 km2), located within the latitudes of 33◦06′N to 43◦01′N and the longitudes of 124◦04′E to 131◦05′E, on a daily basis. The satellite-based image-processing model Surface Energy Balance Algorithms for Land (SEBAL) was adopted and modified to generate the spatial ET data. The SEBAL was calibrated using two years (2012–2013) of measured ETs by an eddy covariance (EC) flux tower at three locations (two in a mixed forest area and one in a rice paddy area). The primary inputs for the model were land surface temperature/emissivity (LST/E), the Normalized Distribution Vegetation Index (NDVI), albedo (Ab) from a Terra Moderate-resolution Imaging Spectroradiometer (MODIS) satellite, a digital elevation model, and wind speed and solar radiation (Rs) from 76 ground-based weather stations. When LST data were unavailable because of clouds and/or snow, the bias-corrected ground temperature measured at the weather stations was used. The NDVI and Ab were used as the monthly average value to maintain relatively stable values rather than using the original time interval data. The determination coefficient (R2) between SEBAL and the flux tower ET was 0.45–0.54 for the two mixed forest towers and 0.79 for the rice paddy tower reflecting the known characteristics of closed and open space ET estimation. The spatial distribution of SEBAL showed that the spatial ET reflected the geographical characteristics, revealing that the ET of lowland areas was higher than that of highland areas.


Introduction
In general, evapotranspiration (ET) is a term used to describe total loss of water from the Earth's surface to the atmosphere by combined processes of evaporation from the open water bodies, bare soil and plant surfaces, etc. and transpiration from vegetation or any other moisture containing living surface [1], and a total loss from ET accounts for approximately 40% out of total amount of water resources in South Korea [2].During the past several decades in South Korea, dramatic climate change has manifested as increased temperature and the number of heavy precipitation.For the spatial and temporal variation of annual precipitation and runoff during 34 years  in Korea, the long-term trends shows the increasing direction in the northern part of Korea and decrease in the southwestern part of Korea [3], although not changed in total amount during a hydrologic year.These changes hinder sustainable water resources management, particularly because of the droughts that arise due to deepening severity and the shortened return period of droughts.The stable securing of water resources by understanding the watershed hydrologic cycle is desired and should be considered as preparation for the future climate change adaptation.
Among the hydrological components, ET is an important component of the watershed hydrological cycle and the most important climatic element, dominating water balance and plant productivity [4].Because ET has a large effect on other related hydrological processes, such as soil moisture dynamics, groundwater discharge, and runoff generation [5], the right understanding in ET variations can provide clues for solving water resources problems.Therefore, the spatial distribution and temporal variation of ET is an important issue in the watershed hydrologic study [6].
Spatial ET measured indirectly from remote sensing (RS) has radically changed the capability of water resources management [7,8].There are various RS-based ET models, such as the Surface Energy Balance Algorithm for Land (SEBAL) [9], the Simplified Surface Energy Balance Index [10], the Surface Energy Balance System [11], mapping evapotranspiration at high resolution with internalized calibration (METRIC) [12], and surface temperature-vegetation index triangle method (Ts-VI triangle model) [13].Among the many approaches, SEBAL has advantages in estimating the surface-energy fluxes from RS data using minimal auxiliary ground-based data and is able to include an automated self-calibration process in each region of interest through the identification of dry and wet pixels and the determination of the near-surface air temperature difference [1,8].
SEBAL can be applied with any satellite sensor having spectral bands, including ASTER (Advanced Spaceborne Thermal Emission and Reflection), CBERS (China-Brazil Earth Resources Satellite), Landsat ETM/TM (Enhanced Thematic Mapper), Terra MODIS (Moderate-resolution Imaging Sepctroradiometer), NOAA-AVHRR (National Oceanic and Atmospheric Administration-Advanced Very High Resolution Radiometer), and RESURS [14].Recently, several studies on SEBAL have used MODIS [8,[15][16][17][18][19][20].The main advantage of MODIS data is the high temporal resolution because such data can be used to estimate energy fluxes at regional, continental and global scales at daily time intervals [21].However, ET estimation of SEBAL is largely subject to the spatial scale effect, which was studied by original analytical models developed by [8].Recently, this effect has been further analyzed by [22], which confirmed that SEBAL is superior to other contexture information based models (e.g., the Ts-VI triangle model) in ET estimation.
In some previous studies, SEBAL turned out that it has several drawbacks: (1) it requires subjective specifications of representative anchor pixels; (2) when SEBAL is applied over mountainous areas, adjustments based on a digital elevation model (DEM) need to be made to account for the lapse rate; and (3) errors in surface temperatures or surface-air temperature differences have great impacts on sensible heat flux estimates [1].Although METRIC based on a similar approach with SEBAL can be used to avoid the limitations of SEBAL [23], METRIC requires computation of alfalfa reference ET for each image and readily available high quality hourly agricultural weather data for ET calculation [24,25].In addition, it was found that SEBAL/METRIC models had high potential for successful ET estimates in the semi-arid US by comparing the derived ET with lysimeter-measured values [24].
In South Korea, SEBAL model was applied in several watershed: Bochung-stream basin [26], Gyeongancheon watershed [27], and Wangsukcheon watershed [28].However, there has not been enough investigation about long-term accuracy of ET estimation in region of temperate climate.In earlier studies [29,30], while SEBAL model was built and validated to estimate daily ET in mixed forest and rice paddy regions, it had some problems including poor efficiency.Therefore, the objective of this study is to estimate daily spatial evapotranspiration using a modified SEBAL algorithm with improved automatic selection for anchor pixel and applied monthly data to enhance the rate of results acquisition using Terra MODIS data for South Korea.The spatial ET results obtained using SEBAL will be validated using two years (2012-2013) of flux tower ET data measured at three locations (two in a forest and one in a rice paddy).

Theory of SEBAL
SEBAL is based on a complete radiation and energy balance along with the resistances for momentum, heat and water vapor transport for each pixel in a given area [9,31].The energy balance equation is explained in 4 parts: where λET is the latent heat flux (W/m 2 ), R n is the net radiation at the surface (W/m 2 ), G is the soil heat flux (W/m 2 ), and H is the sensible heat flux to the air (W/m 2 ).The SEBAL algorithm calculates the net radiation (R n ), soil heat flux (G), and sensible heat flux (H).The latent heat (λET) is obtained via the algorithm described in Figure 1, and the evapotranspiration is computed using the evaporative fraction [λET/ (R n − G)] [32].

Theory of SEBAL
SEBAL is based on a complete radiation and energy balance along with the resistances for momentum, heat and water vapor transport for each pixel in a given area [9,31].The energy balance equation is explained in 4 parts: where is the latent heat flux (W/m 2 ), is the net radiation at the surface (W/m 2 ), is the soil heat flux (W/m 2 ), and is the sensible heat flux to the air (W/m 2 ).The SEBAL algorithm calculates the net radiation ( ), soil heat flux ( ), and sensible heat flux ( ).The latent heat ( ) is obtained via the algorithm described in Figure 1, and the evapotranspiration is computed using the evaporative fraction [ /( − )] [32].When applying the SEBAL model, a crucial point is the selection of two "anchor" pixels (the "hot" and "cold" pixels) over the area of interest.These two pixels are used to determine the differences in temperature between the surface temperature ( ) and air temperature ( ).A linear relationship is assumed between and in the form of where and are the linear relationship constants.To determine these constants, SEBAL uses the two "anchor" pixels for which a value for H can be reliably estimated.is estimated from the land surface temperature (LST) (MOD11_L2) for each pixel; is calculated for either the "hot" or "cold" pixel using the following relationship: When applying the SEBAL model, a crucial point is the selection of two "anchor" pixels (the "hot" and "cold" pixels) over the area of interest.These two pixels are used to determine the differences in temperature between the surface temperature (T s ) and air temperature (dT).A linear relationship is assumed between T s and dT in the form of where a and b are the linear relationship constants.To determine these constants, SEBAL uses the two "anchor" pixels for which a value for H can be reliably estimated.T s is estimated from the land surface temperature (LST) (MOD11_L2) for each pixel; dT is calculated for either the "hot" or "cold" pixel using the following relationship: where H can be calculated for the anchor pixels using meteorological data, ρ is the air density (kg/m 3 ), c p is the specific heat of air (1004 J/kgK), dT is the temperature difference between two heights (K), and r ah is the aerodynamic resistance to heat transport (s/m) for each "cold" and "hot" pixel.The above linear relationship between dT and T s is a major presumption in SEBAL.Several studies indicate that this presumption fits a wide range of conditions [33].Equation ( 3) is developed using dT values for the cold and hot pixels and surface temperature.The cold pixel is used to define the amount of evapotranspiration, through H, occurring from the most vegetated and well-watered areas in the satellite image.A full surface covered crop cultivation or water body is typically used to identify the cold pixels in the area of interest [34].In this study, the automatic selection algorithm of the anchor pixels against the standard is applied as explained in Section 2.2.However, the surface temperature that is used needs to be uniformly adjusted to a common reference elevation for accurate prediction of dT.Otherwise, high elevations that appear to be "cool" may be misinterpreted as having high evaporation.Therefore, SEBAL considers a "lapsed" surface temperature map to compute the surface-to-air temperature differences by assuming that the rate of decrease in surface temperature by the orographic effect is the same as that for a typical air profile.For this reason, DEM data were used for this calculation.The DEM corrected surface temperature is calculated by following equation: where ∆z is the difference of a pixel's elevation from the datum in meters, and is positive if the elevation of a pixel is higher than the datum; and 0.0065 is a constant based on the general statement that air temperature decreases 6.5 • C when elevation increases by 1 km under neutral stability [33].
From the residual in the instantaneous energy-balance equation and the evaporative fraction (Λ), the daily ET (ET 24 , mm/day) was estimated.Λ is notably regular and relatively constant on cloud-free days [35].Thus, its instantaneous value can be taken as the daily mean value, so that the spatial variability in daily ET can be predicted over large scales [20].
where R n,24 is daily net radiation; G 24 is daily soil heat flux; 86,400 is the number of seconds in a 24 h period; and λ is the latent heat of vaporization (J/kg).The latent heat of vaporization allows expression of ET 24 in mm/day.The parameter G 24 can be approximated for vegetative and soil surfaces as zero at the soil surface.This is because, on average, the energy stored in the soil during the daytime is released into the air at night.The latent heat vaporization and R n,24 are defined as [33]: R n,24 = (1 − α) Rs 24 − aτ sw (8) where Rs 24 is 24 h incoming solar radiation; α is albedo; a is a regression coefficient of relationship between net long wave radiation and atmospheric transmissivity at daily scale [36]; and τ sw is the one-way transmittance in clear sky and can be predicted for clear, relatively dry atmospheric conditions using the elevation above sea level in m [37].The coefficient a can be applied differently according to region.In semi-arid conditions, the coefficient a averages 143 with R 2 = 0.80, while the originally suggested value of 110 was found for the Netherlands.In addition, it was applied with 135 in grassland catchments in the Netherlands [38].In this study, 110 was used for the coefficient a and it can be improved in further study to suit this local area.

SEBAL Implementation
In this study, the SEBAL model was applied differently from its typical use [9,33].First, daily input data were required to calculate daily ET, and the input data of some days could be used for estimating results at the time of satellite data acquisition.However, the NDVI and albedo outputs of the Terra MODIS satellite are not available every day, and if we calculated daily NDVI or albedo from other outputs such as surface reflectance, the calculated values would be difficult to use for the model because of noise.Furthermore, albedo products with 16-day temporal granularity are also noisy.The 16-day albedo data were merged into monthly data to redeem fault of the data that had many missing values and outlier values.However, in spite of complementation of the data, the missing values and outlier values were still witnessed.For this reason, the merged albedo data were calibrated [39] according to land use, and the monthly NDVI were considered for consistency.Vegetation gradually changes from dense to thin and a rapid change does not occur in the forest.Thus, monthly albedo and NDVI data could be used for calculating daily ET, and the other data such as LST and meteorological data were obtained as daily data.
Second, the user should search for the hot and cold anchor pixels to set the boundary conditions for pixels exhibiting 0 or nearly 0 ET (hot pixel) and those with maximum ET (cold pixel).In this case, the user should choose the locations of the anchor pixels in each image, but this does not always occur because people can make mistakes as a result of external factors such as psychological and physical states.Furthermore, the criteria for selecting anchor pixels can be changed depending on the user and thus the improper selection of these pixels might cause bias in the sensible heat flux calculation and the resulting ET [40].The METRIC model that departs from the SEBAL model in its use of weather-based reference ET to establish energy balance conditions at a cold pixel attempts to overcome the disadvantage of user dependency through the automatic selection of anchor pixels [41].The METRIC model is based on the same principles and techniques as the SEBAL model [12]; thus, we can also apply the automatic selection system for the SEBAL model.In this study, the automatic selection of anchor pixels in the SEBAL model was performed as follows.For cold pixels, select the top 5% of the highest NDVI, and, among them, select the coldest 15% of T s within an agricultural area of the study area according to land use.For hot pixels, select the lowest 10% of NDVI and then select the hottest 15% of T s within a bare field or an urban area according to land use.The standard for the anchor pixel candidates is referenced by [41], however, the percentage is adjusted slightly from 20% to 15% because there are many pixels in the area of interest and it can affect total run-time of the model.During the period of this study (2012-2013), the maximum value of LST showed the range from about 279.2 K to 321.0 K, while the minimum value of LST indicated the range from around 252.1 K to 295.4 K. Generally, the hot pixels were selected in the range from 277.1 K to 319.1 K, and the cold pixels were included in the range from 253.6 K to 296.4 K in this study.Closer to the summer (from June to August), the temperatures of hot and cold pixels were much higher, and the opposite was shown closer to winter (from December to February).As a result, sub-zero temperatures were selected as cold pixel in some cases in winter season, and the simulated ET tends to be smaller and close to zero in some areas, corresponding to a natural phenomenon.

Study Area and Flux Tower ETs for SEBAL Validation
The study area is the entire area of South Korea, located in northeast Asia, within the latitude and longitude range of 33 • 06 N to 43 • 01 N and 124 • 04 E to 131 • 05 E, respectively.The total area of South Korea is 99,900 km 2 , with more than 65% of the area being forested and 12.2% of the lowland area being crop cultivated.The climate is characterized by four seasons, with heavy precipitation (PCP) during the monsoon summer season (June to August).Under the influence of migratory anticyclones, the spring (March to May) and fall (September to November) seasons are particularly dry.The winter The flux tower sites SMK and CFK (Figure 2) are within Han River basin, the biggest basin in South Korea (total area of 34,406 km 2 ).The SMK tower is located in a steep slope mountainous region within mixed forest of Quercus variaabilis, Quercus mongolica, and Pinuskoraiensis species and with average canopy height of 15 m.The soil has shallow depth with average 50 cm and relatively low retention of soil water contents with monthly average ranging from 11% to 20%.By contrast, the CFK tower is located in a flat farmland region within rice paddy.During the rice growing season from the end of May to the first period of September, the canopy height grows from 5 cm to maximum 1 m with ponding condition.Other periods are not maintained with residue covered surface with relatively high soil water contents (20%-40%).The flux tower sites SMK and CFK (Figure 2) are within Han River basin, the biggest basin in South Korea (total area of 34,406 km 2 ).The SMK tower is located in a steep slope mountainous region within mixed forest of Quercus variaabilis, Quercus mongolica, and Pinuskoraiensis species and with average canopy height of 15 m.The soil has shallow depth with average 50 cm and relatively low retention of soil water contents with monthly average ranging from 11% to 20%.By contrast, the CFK tower is located in a flat farmland region within rice paddy.During the rice growing season from the end of May to the first period of September, the canopy height grows from 5 cm to maximum 1 m with ponding condition.Other periods are not maintained with residue covered surface with relatively high soil water contents (20%-40%).Two flux towers operated by Hydrological Survey Center (HSC) are utilized to monitor the water and carbon cycles between the atmosphere and terrestrial ecosystems based on the EC system.The EC systems for flux monitoring such as water vapor or carbon are mounted 19.2 m high in SMK and 10 m high in CFK.The SMK is an open path system that includes three-dimensional sonic anemometers (CSAT3, Campbell Scientific, Logan, UT, USA; RM-81000, RMYoung, Traverse City, MI, USA) and gas analyzers (LI-7500, LI-COR), while the CFK is a closed path system consisting of three-dimensional sonic anemometers (CSAT3) and gas analyzers (EC-155, Campbell Scientific) [43].
The other site, DMK, operated by K-water (Korea Water Resources Corporation, Daejon, Korea) is located in other basin of Geum River with EL. 688.568 m and 25 m height to avoid the effects of reservoir and canopy height.It is within a mountainous region with mixed forest of soft and hardwood and Larix leptolepis.It has three-dimensional sonic anemometers (CSAT3) and gas analyzers (EC155) installed at 19 m of the tower [44].
The flux towers provide several key meteorological variables as well as flux measurements.The net radiation and other meteorological variables such as air temperature, pressure, humidity and Two flux towers operated by Hydrological Survey Center (HSC) are utilized to monitor the water and carbon cycles between the atmosphere and terrestrial ecosystems based on the EC system.The EC systems for flux monitoring such as water vapor or carbon are mounted 19.2 m high in SMK and 10 m high in CFK.The SMK is an open path system that includes three-dimensional sonic anemometers (CSAT3, Campbell Scientific, Logan, UT, USA; RM-81000, RMYoung, Traverse City, MI, USA) and gas analyzers (LI-7500, LI-COR), while the CFK is a closed path system consisting of three-dimensional sonic anemometers (CSAT3) and gas analyzers (EC-155, Campbell Scientific) [43].
The other site, DMK, operated by K-water (Korea Water Resources Corporation, Daejon, Korea) is located in other basin of Geum River with EL. 688.568 m and 25 m height to avoid the effects of reservoir and canopy height.It is within a mountainous region with mixed forest of soft and hardwood and Larix leptolepis.It has three-dimensional sonic anemometers (CSAT3) and gas analyzers (EC155) installed at 19 m of the tower [44].
The flux towers provide several key meteorological variables as well as flux measurements.The net radiation and other meteorological variables such as air temperature, pressure, humidity and precipitation are measured with net radiometers (CNR2, Campbell Scientific) and automatic weather transmitters (WXT510, Campbell Scientific) at the same measurement heights as the EC systems of each flux tower.Measurements are made at flux towers every 30 min.To reduce the uncertainties associated with flux tower measurements, Quality Assurance (QC) is conducted with respect to planar fit rotation, Webb-Pearman-Leunung (WPL) and spike detection based on statistical approaches [45].In addition, gap-filling methods are applied for missing latent heat flux data using the Mean Diurnal Variation (MDV) or calculating alternative ET based on the Food Agriculture Organization (FAO) 56 Penman-Monteith equation [37,43].
The EC-measured ET at 3 flux towers (SMK, CFK and DMK) was used for the SEBAL validation.Table 1 shows the location and a brief summary of the three flux towers.The EC daily ET from each tower was prepared for 2 years (2012-2013).

SEBAL Input Data
The input data for SEBAL consists of three types of data: RS data, geographic information system (GIS) spatial data, and meteorological data.The RS data were prepared from Terra MODIS satellite products.The Terra satellite, launched in December 1999, views Earth's entire surface every 1 to 2 days.The MODIS sensors on board the Terra satellite acquire data in 36 discrete spectral channels with a spatial resolution of 250 m for visible bands, 500 m for near-infrared bands, and 1000 m for the remaining thermal infrared bands (http://modis.gsfc.nasa.gov/).The data used in this study were retrieved from the Land Processes Distributed Active Archive Center (LP DAAC), including daily LST/emissivity (MOD11_L2) at a spatial resolution of 1000 m, monthly vegetation indices (MOD13A3) at a spatial resolution of 1000 m and 16-day albedo (MCD43A3) at a spatial resolution of 500 m.MOD11_L2 and MOD13A3 data were resampled using majority algorithm in ArcGIS at 500 m spatial resolution to match with the spatial resolution of other data.The majority algorithm is performed to determine the new value of the cell based on the most popular values within the filter window, corresponding 4 by 4 cells in the input space that are closet to the center of the output cell and use the majority of the 4 by 4 neighbors.
The elevation data were rasterized as a DEM with a resolution of 500 m from 1:5000 vector maps supplied by the Korean National Geography Institute (KNGI).The land use data were obtained from the Korean Water resources Management Information System (WAMIS, http://www.wamis.go.kr) with 500 m spatial resolution.SEBAL [6] does not consider land use data.In this study, the land use was used to assist in the selection of the anchor pixels (hot and cold pixels) during the sensible heat flux calculation.Meteorological data from 76 ground weather stations were provided by the Korea Meteorological Administration (KMA, http://www.kma.go.kr) and more information about the station and data is stated in [46].Two categories of meteorological data were prepared during the simulation period: daily ground-measured LST data and daily wind speed observation data.The daily ground measured LST data were observed at 0.0 cm or just slightly under the ground, and daily wind speed measurement was taken in respect to its altitude.
The daily LST was regenerated by gap-filling because the daily MOD11_L2 have missing values because of clouds and atmospheric conditions.The gap-filling was performed each day with bias correction using the following equation recommended by [47,48]: where LST sat is the transformed Terra MODIS satellite LST (MOD11_L2), LST meas is the measured LST of the ground weather station, LST sat is the average satellite LST and LST meas is the average measured LST.

Criteria for the SEBAL Results Validation
Four evaluation criteria were used to evaluate the SEBAL ET with the tower flux ET, i.e., the coefficient of determination (R 2 ), Nash-Sutcliffe efficiency (NSE), Index of Agreement (IOA), Root Mean Square Error (RMSE) and relative RMSE (RMSE%), given by the following equations: where O is the observed value; P is the predicted value; and O and P are the mean observed and mean predicted values, respectively.R 2 ranges from 0 to 1, with higher values indicating a lower error variance; values greater than 0.5 are typically considered acceptable [49,50].NSE ranges between −∞ and 1.0 (1 inclusive), with NSE = 1 being the optimal value.NSE values between 0.0 and 1.0 are generally viewed as acceptable, whereas a value ≤0.0 indicates that the mean observed value is a better predictor than the simulated value, i.e., the performance is unacceptable.IOA ranges from 0.0 to 1.0 and represents the degree to which SEBAL ET is closely related to the tower flux ET.RMSE indicates error in the units (or squared units) of the constituent of interest, which aids in the analysis of the results, with a value of 0 indicating a perfect fit.A detailed description of the evaluation is provided in [51,52].

Correction of the Terra MODIS Data
The Terra MODIS daily LST (MOD11_L2) was regenerated by gap-filling using the 76 ground measured LSTs after bias correction between the MODIS LST and ground LST. Figure 3a shows the comparison of three LSTs of MODIS, ground measured before and after bias correction.Figure 3c Remote Sens. 2016, 8, 983 9 of 20 shows the MODIS LST corrected using Equation ( 9) based on the ground measured LST.For the outlier MODIS LST values in Figure 3a, the values were corrected using the values for the previous and next days.Before correction, the MODIS LST exhibited a deviation of −16.14 K and an R 2 of 0.79 compared with the ground measured LST (Figure 3b).After correction, as shown in Figure 3c, the deviation of MODIS LST decreased by 0.41 K and R 2 increased to 0.99.After bias correction, the gap-filling was conducted to fill the missing values for obtaining the spatial ET distribution.Figure 4 shows the example of MODIS LST before and after the gap-filling of 1 September 2013.First of all, the valid pixel was extracted to point data from satellite LST data in each day (Figure 4b) and subsequently the satellite data were corrected using ground LST of 76 weather stations by bias correction (Figure 4c).One of the reasons why ground LST was considered is the characteristics of climate in the region.Due to the characteristic of monsoon climate, mean cloud cover was higher than 50% generally, and, especially, in summer (from June to August), cloud cover was often observed with more than 90%.Moreover, in some cases, MODIS LST data were clumped on one side.Therefore, to consider spatial distribution of the data, we used ground based LST to fill in missing LSTs, along with bias correction to amend the difference between MODIS and ground LST.Furthermore, to avoid any possible edge effect that can be caused by IDW method, MODIS LST data were extracted from a range of rectangular-shaped area which exceeds the boundary of South Korea.After that, both data were merged (Figure 4d) and interpolated using inverse distance weighted (IDW) method at 500 m spatial resolution (Figure 4e).If MODIS LSTs have few valid data on a particular day, the merged data were replaced by the ground data.Nevertheless, this method is helpful to obtain daily LST data filled with valid data.After bias correction, the gap-filling was conducted to fill the missing values for obtaining the spatial ET distribution.Figure 4 shows the example of MODIS LST before and after the gap-filling of 1 September 2013.First of all, the valid pixel was extracted to point data from satellite LST data in each day (Figure 4b) and subsequently the satellite data were corrected using ground LST of 76 weather stations by bias correction (Figure 4c).One of the reasons why ground LST was considered is the characteristics of climate in the region.Due to the characteristic of monsoon climate, mean cloud cover was higher than 50% generally, and, especially, in summer (from June to August), cloud cover was often observed with more than 90%.Moreover, in some cases, MODIS LST data were clumped on one side.Therefore, to consider spatial distribution of the data, we used ground based LST to fill in missing LSTs, along with bias correction to amend the difference between MODIS and ground LST.Furthermore, to avoid any possible edge effect that can be caused by IDW method, MODIS LST data were extracted from a range of rectangular-shaped area which exceeds the boundary of South Korea.After that, both data were merged (Figure 4d) and interpolated using inverse distance weighted (IDW) method at 500 m spatial resolution (Figure 4e).If MODIS LSTs have few valid data on a particular day, the merged data were replaced by the ground data.Nevertheless, this method is helpful to obtain daily LST data filled with valid data.

Validation of the Energy Balance Component of SEBAL
Figure 5 shows the diurnal variation of the component of energy balance above a well-watered transpiring surface on a cloudless day [37].Another study [53] reported that the net radiation can range from 100 to 700 W/m 2 and provided the / ratio at various surfaces and, in this study, the component of energy balance was validated as the proper ratio proportion.According to the precedence study [23], the net radiation values were estimated to be in the range of 400 to 800 W/m 2 , however, the values tended to be slightly overestimated approximately 100 W/m 2 .The other components, soil heat flux and sensible heat flux, showed a suitable range compared with Figure 5.The soil heat flux is in the range from 0 to 100 W/m 2 at the three flux towers, and the sensible heat flux ranges from approximately 50 to 400 W/m 2 at SMK and CFK and from 100 to 500 W/m 2 at DMK.The difference in the sensible heat flux of DMK compared to those of the other towers is assumed to be due to the effect of elevation.

Validation of the Energy Balance Component of SEBAL
Figure 5 shows the diurnal variation of the component of energy balance above a well-watered transpiring surface on a cloudless day [37].Another study [53] reported that the net radiation can range from 100 to 700 W/m 2 and provided the G/R n ratio at various surfaces and, in this study, the component of energy balance was validated as the proper ratio proportion.According to the precedence study [23], the net radiation values were estimated to be in the range of 400 to 800 W/m 2 , however, the values tended to be slightly overestimated approximately 100 W/m 2 .The other components, soil heat flux and sensible heat flux, showed a suitable range compared with Figure 5.The soil heat flux is in the range from 0 to 100 W/m 2 at the three flux towers, and the sensible heat flux ranges from approximately 50 to 400 W/m 2 at SMK and CFK and from 100 to 500 W/m 2 at DMK.The difference in the sensible heat flux of DMK compared to those of the other towers is assumed to be due to the effect of elevation.

Valibration of SEBAL
As mentioned above, one of the SEBAL algorithms estimates using the instantaneous ET.In the SEBAL model, was calculated from , and .Equation ( 8) calculates , under clear sky conditions using .However, if the day of the satellite image is known to have had some cloudiness during periods preceding or following the time of the image, then one should use a local (ground-based) measured value for 24-h solar radiation ( ) in place of [33].During the

Valibration of SEBAL
As mentioned above, one of the SEBAL algorithms estimates ET 24 using the instantaneous ET.In the SEBAL model, ET 24 was calculated from R n,24 and G 24 .Equation ( 8) calculates R n,24 under clear sky conditions using Ra 24 .However, if the day of the satellite image is known to have had some cloudiness during periods preceding or following the time of the image, then one should use a local (ground-based) measured value for 24-h solar radiation (Rs) in place of Ra 24 τ sw [33].During the simulation period, mean annual cloud cover was higher than 50% (53% and 50% in 2012 and 2013, respectively) as a result of the monsoon climate.For this reason, we find that Ra 24 was unsuitable for South Korea and instead used Rs.The SEBAL ET was estimated using both Ra 24 and Rs, and the scatter plots are shown in Figure 6.Although the SEBAL ET values of Ra 24 were overestimated compared to the flux tower ET at all three sites, the ET of the Rs exhibited similar results because the Rs data incorporated the cloud cover through observation from ground weather stations.
the scatter plots are shown in Figure 6.Although the SEBAL ET values of were overestimated compared to the flux tower ET at all three sites, the ET of the exhibited similar results because the data incorporated the cloud cover through observation from ground weather stations.Table 2 shows the statistical summary of the SEBAL model results.The ET from Ra 24 data indicates that the results were unsatisfactory based on the efficiency criteria.The R 2 between the flux tower and the ET of Ra 24 data during 2012 was 0.26, 0.52, and 0.22 at SMK, CFK, and DMK, respectively.The other statistical results consisted of NSEs of −4.91 (SMK), −0.34 (CFK), and −1.54 (DMK); IOAs of 0.48 (SMK), 0.77 (CFK), and 0.58 (DMK); and RMSEs of 2.10 mm/day (SMK), 1.46 mm/day (CFK), and 1.03 (DMK).In contrast, the ET from Rs data exhibited a reasonably good fit between the flux tower and simulated ET.The R 2 value was 0.53, 0.77, and 0.53 at SMK, CFK, and DMK, respectively.The other statistical results consisted of NSEs of 0.41 (SMK), 0.67 (CFK), and 0.41 (DMK); IOAs of 0.82 (SMK), 0.92 (CFK), and 0.84 (DMK); and RMSEs of 0.67 mm/day (SMK), 0.73 mm/day (CFK), and 0.53 (DMK).Both the ET using Ra 24 data and Rs data indicated that the CFK site produced the most reasonable results.

Application of SEBAL
The SEBAL model was used to estimate the ET from 2012 to 2013. Figure 7 shows the daily time series and scatter plots for the SEBAL ET and the flux tower ETs for two years, and Table 3 shows the statistical summary of the SEBAL ET and the flux tower ETs at three locations.Two flux towers, SMK and DMK, are located in mixed forest areas (closed), whereas CFK is located in a rice paddy area (open).The total amount of ET was found to vary depending on the land use type.In the rice paddy area (CFK), the total flux tower ET was 496.1 mm in 2012 and 467.8 mm in 2013, whereas in the forest areas (SMK and DMK), the flux tower ET was approximately half as large, with values of 311.7 mm and 243.5 in 2012 and 293.8 mm and 255.8 mm in 2013 for SMK and DMK, respectively.The SEBAL ET showed similar patterns as the flux tower ET: the ET of CFK is approximately twice as high as those of SMK and DMK.The maximum ET at SMK during the simulation period was 3.7 and 3.4 mm/day based on the flux tower and SEBAL, respectively; the maximum ET at DMK was 3.6 mm/day and 3.7 mm/day based on the flux tower and SEBAL, respectively; and the maximum ET at CFK was 5.2 and 5.3 mm/day based on the flux tower and SEBAL, respectively.The SEBAL ET was underestimated by 119.5, 110.5, and 45.3 mm/year at SMK, CFK, and DMK, respectively, for the two-year period.The average R 2 values for SMK, CFK, and DMK were 0.53, 0.78, and 0.59, respectively.The NSEs and IOAs for SMK, CFK, and DMK were 0.34, 0.71, and 0.46 and 0.81, 0.92, and 0.86, respectively.Moreover, the RMSE values for SMK, CFK, and DMK were 0.65, 0.66, and 0.58 mm/day, respectively.Although average R 2 values for SMK, CFK, and DMK were greater than 0.5, which were considered acceptable [52], the results could be improved even further by optimizing coefficients including the coefficient a in Equation (8).

Comparison of the Monthly Data at the Three Flux Tower Locations
This study was conducted to calibrate the satellite data for estimating daily ET. Figure 8 shows the monthly data of NDVI, albedo, the flux tower ET, and SEBAL ET at the three flux tower locations.In South Korea, rice is planted from the middle ten days of May to June and harvested from August to September.Whereas the NDVI of the rice paddy area demonstrates the growth and development pattern of the plant from May to September, the NDVI of mixed forest areas exhibits a different pattern, namely, the value higher than 0.5 from April to November.In other words, deciduous trees in the forest have leaves in April and then shed their leaves in September.In the winter season (particularly December), the NDVI of DMK tends to be lower than that of CFK.This result appears to have been influenced by the elevation above sea level.DMK is located at an elevation of 689 m, and the temperature is colder than that at CFK, which is located at an elevation of 141 m.It is assumed that weather conditions caused by the low temperature have a negative effect on vegetation vitality, such as the effects of snow and frost, thereby resulting in low NDVI values.In the summer of 2013, the NDVI of SMK was lower than the value during 2012 at the same site.The NDVI is thought to be affected by rainfall (Figure 8a).In all three sites, the albedo exhibits rising and falling patterns according to the season, similar to the NDVI.In the winter season (December to February), the albedo shows rapidly rising patterns in all sites, and it is assumed that snow has an effect on the albedo

Comparison of the Monthly Data at the Three Flux Tower Locations
This study was conducted to calibrate the satellite data for estimating daily ET. Figure 8 shows the monthly data of NDVI, albedo, the flux tower ET, and SEBAL ET at the three flux tower locations.In South Korea, rice is planted from the middle ten days of May to June and harvested from August to September.Whereas the NDVI of the rice paddy area demonstrates the growth and development pattern of the plant from May to September, the NDVI of mixed forest areas exhibits a different pattern, namely, the value higher than 0.5 from April to November.In other words, deciduous trees in the forest have leaves in April and then shed their leaves in September.In the winter season (particularly December), the NDVI of DMK tends to be lower than that of CFK.This result appears to have been influenced by the elevation above sea level.DMK is located at an elevation of 689 m, and the temperature is colder than that at CFK, which is located at an elevation of 141 m.It is assumed that weather conditions caused by the low temperature have a negative effect on vegetation vitality, such as the effects of snow and frost, thereby resulting in low NDVI values.In the summer of 2013, the NDVI of SMK was lower than the value during 2012 at the same site.The NDVI is thought to be affected by rainfall (Figure 8a).In all three sites, the albedo exhibits rising and falling patterns according to the season, similar to the NDVI.In the winter season (December to February), the albedo shows rapidly rising patterns in all sites, and it is assumed that snow has an effect on the albedo (Figure 8b).Oke [54] and Ahrens [55] reported that the general albedo value of various surfaces and glacier ice are in the range of 0.2 to 0.4.
The monthly ET ranged from 0.9 to 92.6 mm at the flux towers and 0.0 to 90.7 mm from SEBAL (Figure 8c,d).For both sets of data, the ET of CFK was higher than those of the other sites and reached a maximum value in June and then decreased in July and September.In August, ET increased again for both of the datasets in 2013 because of the rainy season in South Korea.At SMK and DMK, the SEBAL ET had similar ranges of 0.0 to 44.0 mm and 0.0 to 43.8 mm, respectively.However, the flux towers at SMK and DMK exhibited different ET patterns of 1.9 to 57.7 mm and 0.9 to 63.0 mm, respectively.The flux tower ET of DMK was missing in January and February.The DMK had a smaller range of ET than those the other sites; thus, we consider that this flux tower was affected by surface elevation.Moreover, the monthly SEBAL ET indicated that the value was changed by the variability of the NDVI and albedo.When the NDVI increased with the change in season from winter to summer, the SEBAL ET increased with an analogous pattern.Meanwhile, the SEBAL ET decreased when the albedo increased in winter, exhibiting an inverse correlation.
Figure 9 shows the spatial distribution of annual SEBAL ET in 2012 and 2013, respectively.The total amount of ET showed similar values in both years.However, the spatial distributions showed distinct difference according to the locations.In Han River watershed (Figure 9), the annual ET showed higher values in 2012 than in 2013.This was because the area in 2013 suffered severe drought, and the annual precipitation was about 35%~50% compared to the long term average from 1973 to 2015 [56].Except the 2012 Han River watershed, the annual ET of other areas such as Guem River and Nakdong River watersheds showed higher values in 2013 than in 2012.Looking at Figure 9, the amounts of ET are different according to land use.In the big city areas of Seoul, Daejeon, Daegu, Gwangju, and Busan, the ET was smaller than the other areas.In addition, the spatial distribution of SEBAL showed the spatial ET reflected the geographical characteristics, revealing that the ET of lowland areas was higher than the ET of highland areas.The monthly ET ranged from 0.9 to 92.6 mm at the flux towers and 0.0 to 90.7 mm from SEBAL (Figure 8c,d).For both sets of data, the ET of CFK was higher than those of the other sites and reached a maximum value in June and then decreased in July and September.In August, ET increased again for both of the datasets in 2013 because of the rainy season in South Korea.At SMK and DMK, the SEBAL ET had similar ranges of 0.0 to 44.0 mm and 0.0 to 43.8 mm, respectively.However, the flux towers at SMK and DMK exhibited different ET patterns of 1.9 to 57.7 mm and 0.9 to 63.0 mm, respectively.The flux tower ET of DMK was missing in January and February.The DMK had a smaller range of ET than those the other sites; thus, we consider that this flux tower was affected by surface elevation.Moreover, the monthly SEBAL ET indicated that the value was changed by the variability of the NDVI and albedo.When the NDVI increased with the change in season from winter to summer, the SEBAL ET increased with an analogous pattern.Meanwhile, the SEBAL ET decreased when the albedo increased in winter, exhibiting an inverse correlation.
Figure 9 shows the spatial distribution of annual SEBAL ET in 2012 and 2013, respectively.The total amount of ET showed similar values in both years.However, the spatial distributions showed distinct difference according to the locations.In Han River watershed (Figure 9), the annual ET showed higher values in 2012 than in 2013.This was because the area in 2013 suffered severe drought, and the annual precipitation was about 35%~50% compared to the long term average from 1973 to 2015 [56].Except the 2012 Han River watershed, the annual ET of other areas such as Guem River and Nakdong River watersheds showed higher values in 2013 than in 2012.Looking at Figure 9, the amounts of ET are different according to land use.In the big city areas of Seoul, Daejeon, Daegu, Gwangju, and Busan, the ET was smaller than the other areas.In addition, the spatial distribution of SEBAL showed the spatial ET reflected the geographical characteristics, revealing that the ET of lowland areas was higher than the ET of highland areas.

Conclusions
This study estimated the spatial ET of South Korea during 2012 and 2013 using the SEBAL model and MODIS satellite data; the results were evaluated by flux tower measurement data at three sites (SMK, CFK and DMK).Daily MODIS data, MOD11_L2, were used after the gap-filling process to mitigate the missing data using ground station measured LST data.The other satellite data were the monthly albedo and NDVI at a spatial resolution of 500 m.The SEBAL results reflected the geographical characteristics, namely, the ET of lowland areas was higher than that of highland areas.The primary results are summarized as follows:

Conclusions
This study estimated the spatial ET of South Korea during 2012 and 2013 using the SEBAL model and MODIS satellite data; the results were evaluated by flux tower measurement data at three sites (SMK, CFK and DMK).Daily MODIS data, MOD11_L2, were used after the gap-filling process to mitigate the missing data using ground station measured LST data.The other satellite data were the monthly albedo and NDVI at a spatial resolution of 500 m.The SEBAL results reflected the geographical characteristics, namely, the ET of lowland areas was higher than that of highland areas.The primary results are summarized as follows: (1) To fulfill gap-filling of MODIS LST, 76 ground measured LSTs were used for bias correction, and the R 2 increased by 0.20 after correction.After that, both data were merged and interpolated using IDW method.If MODIS LSTs have few valid data on a particular day, the merged data were replaced using ground data.On some days, the model calculated zero ETs.There are still many things to be improved.For example, the coefficient a with 110 in Equation ( 8) should be calibrated with more data of measured R n,24 .(6) During winter (particularly December), the NDVI of DMK tended to be lower than that of CFK.This result is apparently influenced by the elevation above sea level.DMK is located at an elevation of 689 m, and the temperature is colder than that at CFK, which is located at an elevation of 141 m.It is assumed that weather conditions caused by low temperature have a negative effect on vegetation vitality, such as the effects of snow and frost, leading to low NDVI values.
This study showed some possibility to apply SEBAL for South Korea spatial ET information by reflecting meteorological and geological conditions.However, the application periods (2012 and 2013) should first be increased, and then validated with flux tower data.In addition, better mining of satellite data and the SEBAL equation constants are necessary to study with observed data.In future study, our method can be improved by developing a regression model between the LST measured at the station and MODIS LST during available days and then this relationship can be used to predict LST at the stations during missing days to fill in gaps in LST.

Figure 1 .
Figure 1.Flowchart of the Surface Energy Balance Algorithm for Land (SEBAL) model.

Figure 1 .
Figure 1.Flowchart of the Surface Energy Balance Algorithm for Land (SEBAL) model.

20 (
December to February) season has a frigid climate compared with other climates at the same latitude due to the influence of the Siberian air mass.The mean annual temperature ranges between 10 • C and 15 • C, with January and August being the coldest (from −6 • C to 3 • C) and hottest (from 23 • C to 26 • C) months, respectively.The mean annual PCP in South Korea is approximately 1300 mm (60%-70% of the annual average PCP rate during the monsoon season) [42].

Figure 2 .
Figure 2. Location of the weather stations and flux towers in South Korea.

Figure 2 .
Figure 2. Location of the weather stations and flux towers in South Korea.

Figure 3 .
Figure 3.Time series and scatter plots of the MODIS and ground measured LST in 2012: (a) comparison of the MODIS and ground measured LSTs; (b) scatter plot of the data before calibration; and (c) scatter plot of the data after calibration.

Figure 3 .
Figure 3.Time series and scatter plots of the MODIS and ground measured LST in 2012: (a) comparison of the MODIS and ground measured LSTs; (b) scatter plot of the data before calibration; and (c) scatter plot of the data after calibration.

Figure 4 .
Figure 4. Filled LST data production process at 1 September 2013: (a) original LST data; (b) valid pixel data extracted from the MODIS LST data; (c) ground measured LST; (d) merged LST data between the MODIS LST and ground measured LST; and (e) gap-filled LST data.

Figure 4 .
Figure 4. Filled LST data production process at 1 September 2013: (a) original LST data; (b) valid pixel data extracted from the MODIS LST data; (c) ground measured LST; (d) merged LST data between the MODIS LST and ground measured LST; and (e) gap-filled LST data.

Figure 5 .
Figure 5.Time series of energy balance component at the three flux towers: (a) SMK; (b) CFK; and (c) DMK.

Figure 5 .
Figure 5.Time series of balance component at the three flux towers: (a) SMK; (b) CFK; and (c) DMK.

Figure 6 .
Figure 6.Scatter plots (January to December 2012) of SEBAL simulation ET using extraterrestrial radiation (left) and solar radiation (right) with flux towers at: (a) SMK; (b) CFK; and (c) DMK.Table 2 shows the statistical summary of the SEBAL model results.The ET from data indicates that the results were unsatisfactory based on the efficiency criteria.The R 2 between the flux tower and the ET of data during 2012 was 0.26, 0.52, and 0.22 at SMK, CFK, and DMK, respectively.The other statistical results consisted of NSEs of −4.91 (SMK), −0.34 (CFK), and −1.54 (DMK); IOAs of 0.48 (SMK), 0.77 (CFK), and 0.58 (DMK); and RMSEs of 2.10 mm/day (SMK), 1.46 mm/day (CFK), and 1.03 (DMK).In contrast, the ET from data exhibited a reasonably good fit

( 2 )
In the time series of the energy balance equation component, the net radiation was overestimated by approximately 100 W/m 2 compared to the results of preceding papers in the summer season.The other components, soil heat flux and sensible heat flux, were within a valid range, although the sensible heat flux of DMK is in a different range because of the effect of elevation.(3)When we applied solar radiation (Rs) instead of Ra 24 τ sw , the R 2 of SEBAL model was improved by almost double in forest area, and also increased in rice paddy area from 0.52 to 0.77.This approach may be useful with almost cloudy weather condition except short spring and autumn periods of clear sky in a year.(4) The model results varied depending on land use type.In CFK, a rice paddy area, the total flux tower ET was 496.1 mm in 2012 and 467.8 mm in 2013, whereas, in SMK and DMK, which are both forest areas, the flux tower ET was approximately half of the amount at 311.7 and 243.5 mm, respectively, in 2012 and 293.8 and 255.8 mm, respectively, in 2013.SEBAL ET exhibited similar patterns with the flux tower ET: the ET of CFK is approximately twice as high as those of SMK and DMK.(5) The SEBAL results showed average R 2 values of 0.53 (SMK), 0.78 (CFK), and 0.59 (DMK), and IOAs of 0.81 (SMK), 0.92 (CFK), and 0.86 (DMK), although we used some monthly data in the model application.The LST was the key factor besides satellite data to estimate daily temporal variability.

Table 1 .
Summary of the three flux towers.

Table 2 .
Comparison SEBAL results of Ra 24 and Rs ET at three flux tower locations in 2012.

Table 3 .
Statistical summary of the SEBAL ET and flux tower ET for two years (2012 to 2013).