Soil moisture derived using two apparent thermal inertia functions over Canterbury, New Zealand

Abstract The near-surface soil moisture (SM) is an important property of the soil that can be studied from satellite remote sensing observations over a large spatial domain. This research provides an estimate on the accuracy of SM retrieved from satellite land surface temperature (LST) observations over the Canterbury Plains, New Zealand. The apparent thermal inertia (ATI) method with two approaches (ATI1 and ATI2) was applied to derive the near-surface SM from the moderate resolution imaging spectroradiometer (MODIS) LST product. The in-situ measurements of SM and rainfall data at six sites across the study area were used as reference. The analysis was conducted over two periods, a short period of four months and a longer period of three years. SM simulations by the weather research and forecasting (WRF) model were used in the analysis for the shorter period. Overall, SM based on ATI2 showed a slightly higher correlation with the in-situ measurements ( ρ ¯ = 0.66 ) than ATI1 ( ρ ¯ = 0.63 ). The correlation, in general, was higher for the WRF simulations ( ρ ¯ = 0.81 ). Both functions performed better during summer compared to winter, but overall, ATI2 showed lower mean errors ( ME ≈ − 15   m 3 · m − 3 volumetric SM) compared to ATI1 ( ME ≈ − 20   m 3 · m − 3 ) at most of the sites. Additionally, seasonal variations of SM were better detected by ATI2 than ATI1, and the effects of precipitation were detected on more occasions by the ATI2 function. We conclude that ATI2 function can be used to estimate the near-surface SM over a large area from the MODIS LST time series if a few representative reference stations are available.


Introduction
The near-surface soil moisture (SM) affects diurnal change of surface temperature, and is a key variable in computing several parameters of the land energy and water budget. 1 Observations of spatially distributed SM are essential for a large range of hydrological, climate, and agricultural applications. 2 Satellite observations in the visible and near-infrared (VNIR), thermal infrared (TIR), and microwave (MW) regions of the electromagnetic radiation can be used to derive the near-surface SM 3,4 with a varying depth depending on the surface type and the energy employed by the sensor. Land surface temperature (LST) is the parameter measured by TIR observations. One of the algorithms based on TIR datasets is the thermal inertia (TI) method. 5 The apparent thermal inertia (ATI), a quantification of the effect of soil TI on soil surface temperature, is solely based on remotely sensed observations. 3,4,6 ATI makes use of the difference between day and night LST as well as surface albedo information derived from VNIR observations. 7 Methods for retrieving the near-surface SM using ATI mainly differ in their approach to estimate the diurnal cycle of LST. Ref. 8 developed a model (known as XC) to derive the real TI from the advanced very high resolution radiometer (AVHRR) day (14:37 local solar time) and night (04:44) observations and found that TI is directly proportional to ATI. To compute the real TI, the XC method required an extra parameter measured on the ground, which was the time of the maximum temperature in daytime obtained from a meteorological station. Ref. 9 developed a method based on three AVHRR LST observations per day (02:30, 7:30, and 14:30 local solar time) to retrieve TI without the need for a ground measured parameter. The same authors developed another model, 10,11 named FTA (four temperatures algorithm), using four AVHRR thermal observations (02: 30, 7:30, 14:30, and 19:30) and compared their method with XC. They found that FTA is better in presenting the diurnal cycle of soil surface temperature as it uses more LST observations over a 24-h period. Ref. 6 also used two LST observations (day and night) from Meteosat, but to estimate ATI rather than TI. They approximated the topsoil saturation index from ATI and scaled it using the maximum and minimum long-term SM measurements on the ground to obtain the near-surface SM content (we will call this ATI1 function). Ref. 4 modified FTA to be applied on the moderate resolution imaging spectroradiometer (MODIS) LST observations from Terra and Aqua day and night (∼1:30, 10:30, 13:30, and 22:30) overpasses (we will call it ATI2 function).
The LST product is one of many datasets derived from day and night observations of MODIS onboard Terra and Aqua satellites for more than a decade. MODIS BRDF/Albedo product is another dataset, which is useful as a necessary input for ATI calculations. As a property of land surface, SM can also be estimated based on surface moisture fluxes; therefore, coupled land atmospheric models with regional scales are also used to simulate SM across spatial extent. The weather research and forecasting (WRF) community model is a mesoscale regional model, which can be used to simulate SM over scales ranging from meters to thousands of kilometers 12 when it is coupled with one of the available land surface schemes in the model. Using numerically modeled SM in combination with a remote sensing approach can be helpful to fill the gaps in the latter due to cloud cover. It also provides the opportunity for intercomparison of the two outputs over spatial domain. This is critically important on areas where ground truth data are not available (such as rugged alpine or densely forested parts of the South Island of New Zealand).
The objective of this research, therefore, is to assess the relationship between SM derived from the MODIS LST using the ATI method, simulated by the WRF model and measured on the ground, and to compare the performance of ATI1 and ATI2 functions. The study is aimed to find out if a more detailed representation of the LST diurnal cycle in the ATI model leads to an improvement in the near-surface SM retrievals. Based on our review of the literature, there have been limited attempts previously to compare different ATI methods and to evaluate their performance in relation to rainfall.

Study Area
The study area covers part of the Canterbury Plains in the South Island of New Zealand (Fig. 1), approximately centered at 43°37′ S and 172°11′ E. Surface elevation of the study area ranges from a few meters near the coast to about 200 m above sea level (asl) near the Southern Alps in the West. The Canterbury Plain is the largest area of flat land in New Zealand dominated by agricultural and farming land use. Widespread use of irrigation is an ever increasing need in the region to sustain productivity. For the same reason, study of soil water content over a large spatial extent is a prime concern in this region. Only satellite data and modelling approaches can provide an adequate coverage of SM variability over this extensive spatial domain. To obtain consistent results, however, remotely sensed data and modeled parameters need to be validated based on ground measurements.
Measurements of SM at six sites across the study area ( Fig. 1) were used as reference. The land cover (LC) types of these sites were dominantly grass mixed with trees, irrigated crops, urban developed areas, and barren/fallow land ( Fig. 2 and Table 1).

Data
Data used in this paper include the in-situ measurements of SM and rainfall, the remotely sensed LST and albedo products, and the WRF simulations.

In-situ Measurements
The in-situ measurements of SM (20-cm depth) over a period of three years (January 01, 2010, to December 31, 2012) recorded at six automatic stations across the study area ( Fig. 1) are used in this research. These data were downloaded from the national climate database-Cliflo (http:// cliflo.niwa.co.nz) of New Zealand, which is maintained by the National Institute of Water  and Atmospheric Research (NIWA, http://www.niwa.co.nz). These data were recorded with an hourly rate, but were downgraded to a daily rate to match the daily frequency of ATI SM derivations. Daily rainfall data, accumulated over a period of 24 h, were also downloaded from these stations for the same period to be used in the analysis. Additional in-situ measurements of the near-surface SM (2 to 5 cm depth) over five LC types were also collected in the field; however, these measurements turned out to be highly affected by solar radiation; therefore, they were not used in the analysis.

MODIS LST Product
This product is an operational dataset derived from MODIS thermal observations onboard Terra and Aqua satellites. These observations are collected during day and night overpasses of Terra and Aqua (approximately at 1:30, 10:30, 13:30, and 22:30) with different viewing angles. The LST product is derived from bands 31 and 32, which are observed at 10.78-11.28 μm and 11.77-12.27 μm spectral ranges, respectively. Theoretical background and technical details of the algorithms and the procedure for extraction of LST from MODIS thermal bands are beyond the scope of this paper, but can be accessed in the literature. [13][14][15]

MODIS Combined BRDF/Albedo Product
The combined albedo product (MCD43B3) is generated from a combination of MODIS-Terra and MODIS-Aqua albedo products (MOD43B3 & MYD43B3) accumulated over a 16-day period. MOD43B3 and MYD43B3 are Level 3 global products with 1-km resolution mapped into a sinusoidal grid. As a result, MCD43B3 combined albedo product is an 8-daily dataset with the same resolution and geometric specifications 16,17 as the two preceding products. This dataset contains two types of albedos: a "black-sky albedo" (BSA), which is the directional hemispherical reflectance that integrates the BRDF over the exitance hemisphere for a single irradiance direction, and a "white-sky albedo" (WSA) which is the bi-hemispherical reflectance that integrates the BRDF over all viewing and irradiance directions. 18 Data for each of the BSA and WSA categories are derived from MODIS bands 1-7. Three more integrated broadband albedos from the visible (0.3-0.7 μm), near-infrared (0.7-3.0 μm), and shortwave (0.3-5.0 μm) regions 17,19 for each of the BSA and WSA categories are also provided in this dataset. Broadband albedo is the ratio of radiant energy scattered upward and away from the surface in all directions to the downwelling irradiance incident upon the surface 18 from all directions. The actual albedo on the ground, or "blue-sky albedo" can be estimated as a sum of BSA and WSA weighted by the proportions of direct and diffuse solar radiations arriving at the ground. 19 We interpolated both BSA and WSA fields in the 8-daily combined albedo products to daily values using a linear interpolation function. In the next step, we calculated the blue-sky albedo based on the sum of the shortwave broadband BSA and WSA weighted by 0.34 and 0.66, respectively. The advantage of this approach, rather than using a constant value for albedo in ATI calculations (which was practiced by other references 8,10 ), is taking the seasonal variability of the surface albedo into account.

WRF Simulations
Surface layer SM simulations from the WRF model for a period of four months (July 31 to December 01, 2010) were used in this paper. We used the Noah land surface model (LSM) as a land parameterization scheme. 20 Noah LSM is one of the multiple land surface schemes coupled with the WRF modeling system. 21 The initial conditions for SM simulations in the model are defined by the National Center for Environmental Prediction (NCEP) global final analysis (FNL) data with a spatial resolution of 1 × 1 deg (∼80 × 111 km in our study area) and 6-h frequency. 12 The input variables from the NCEP FNL data include (but are not limited to) soil moisture/water content, soil temperature, precipitation, heat flux, humidity, surface winds, and land cover (see http://rda.ucar.edu/datasets/ds083.2). The Noah LSM only provides surface heat and moisture fluxes as lower boundary conditions to the coupled atmospheric (i.e., the WRF) model. 22 These fluxes are then transported throughout the boundary layer, and interact with other model physics involving cloud, radiation, and precipitation processes. 22 In the Noah LSM, the green vegetation fraction (Fg) parameter, which is defined as the coverage of vegetation over an area unit or a pixel, 12 is used to differentiate moisture flux from vegetation. As a result, the WRF model coupled with the Noah LSM is able to separate soil and vegetation moisture fluxes and also accounts for precipitation effects in SM approximation. To calculate SM, a prognostic diffusion equation for the volumetric SM content is used. 23 Sensitivity of the WRF model to different land surface parameterization schemes, such as the Noah LSM, is discussed in the literature, 12,21 which includes SM simulations. Noah LSM simulates the near-surface SM in four predefined layers (with a top-down sequence of 0.05, 0.25, 0.70, and 1.5 m depth). Outputs for the top most layer (5-cm depth) were used in our analysis. The first 24 h of the simulations were discarded from the analysis as the spin-up period required by the model to reach a balanced state with the boundary conditions. Grid spacing of the simulations was set to 5 km for the first and 1 km for the second domain of the model to match the spatial resolution of the MODIS LST. The time interval of the model outputs was set to a 30-min rate for all simulations.

Thermal Inertia Approach for Soil Moisture Estimation
TI, defined as the resistance of a material to change in temperature, 24 is calculated based on the knowledge of thermal conductivity and density of the near-surface soil layer Eq. (1) is the soil bulk density, and C (J · kg −1 · K −1 ) is the soil heat capacity. 3,25,26 Water bodies have a higher TI than dry soils and rocks, and exhibit a lower diurnal temperature amplitude (DTA), therefore, when soil water content increases, DTA of the soil will decrease. 6 This property can be exploited to derive the amount of moisture content in the upper soil from remotely sensed LST. It must be noted that in our investigation area, TI is also affected by the low to medium (mainly grass) vegetation ( Fig. 1). The ATI method based on two or more daily LST observations and one daily albedo data, along with a priori knowledge about the acquisition date and the geographic latitude, is used to derive the near-surface SM from remotely sensed observations Eq. (2). ATI is defined as 3,4,6,24 where α 0 is the broadband albedo, DTA is derived from two or more daily LST observations, and S is the solar correction factor defined as where ϑ is the local latitude and φ is the solar declination, which is calculated for each day of the year. 4 If only two observations representing coolest and warmest LST values over a 24-h period are used, the DTA will be simply the difference between the diurnal warm and cool temperatures (ΔLST), defined as 6,24 where LST day and LST night are daily LST observations at peak warm and cool hours, respectively. LST observations from MODIS-Aqua with approximate overpass times at 1:30 am and pm (representing LST night and LST day ) have been often used in the literature to calculate DTA. 6 MODIS-Terra observations with approximate overpass times at 10:30 am and pm (representing LST day and LST night , respectively) were also used by others. 24 A second method to calculate DTA is to apply four LST values from a combination of daily MODIS-Terra and MODIS-Aqua observations 4 In Eq. (5), n ¼ 4 for four MODIS daily observations, ω ¼ 2π∕ð24 Ã 60 Ã 60Þ rad s −1 is the angular velocity of the Earth's rotation, t i [i ¼ 1; : : : ; 4] is the overpass time, T i is the surface temperature (i.e., LST) at each overpass time, and ψ is the phase angle.
The near-surface SM can be then calculated from the ATI values for every pixel 6 after checking for possible outliers. After calculation of ATI, Ref. 6 used the minimum and maximum ATI values over each pixel over time (t) to calculate the remotely sensed topsoil saturation index (SMSI RS ) using Eq. (6): where ATI ðtÞ is the ATI value at time t, ATI min is the minimum, and ATI max is the maximum ATI values in the entire time series of the pixel under analysis. Finally, Ref. 6 derived volumetric soil moisture content (SMC) based on the SMSI RS values where t is time, SM res is the residual volumetric SMC, and SM sat is the saturation volumetric SMC. Since SM res and SM sat can only be determined under laboratory conditions, field measurements are unlikely to be equal to extremely high or low laboratory values. 6 For the same reason, SM res and SM sat in practice are substituted by the minimum (SM min ) and maximum (SM max ) in-situ measurements of SM over the entire period, respectively (Table 1).
In this paper, ATI as given in Eq. (2) was calculated with two different approaches: ATI1 using Eq. (4) and ATI2 using Eq. (5), each of which takes a different solution to estimate DTA. In the next step, Eq. (6) is used to derive SMSI RS at each site. Finally, we scaled SMSI RS values to the respective SM min and SM max measurements at each site using Eq. (7). Table 1 Geographic coordinates of the in-situ soil moisture (SM) and rainfall measurement stations alongside the minimum (SM min ) and maximum (SM max ) measured SM at each site during the measurement period (3 years: January 2010-December 2012), and the dominant land-cover around each site. The analysis was conducted in two temporal scales: short term (four months: August 01 to December 01, 2010) and long term (three years: January 01, 2010, to December 31, 2012). The first time scale was chosen to identify the effects of rainfall on ATI SM retrievals, as well as the WRF SM simulations, and to evaluate the difference in seasonal trends in comparison to the in-situ measurements. Also, since model simulations of SM with a 1-km spatial resolution are time consuming and memory intensive (especially when the time period is longer than a few months), the analysis involving the WRF model was conducted only for the short-term period. The long-term analysis was necessary to evaluate the quality of ATI SM retrievals in different seasons and to find out which ATI function is more able to detect seasonal variability of SM.

Statistical Methods
The Euclidean distance (d E ) was applied to compare the difference between the time series of SM from ATI functions with the in-situ measurements Eq. (8) where ts 1 and ts 2 are two time series to be compared (e.g., ATI1 and in-situ SM series), and n is the total number of coincident time points in both series. If ts 1 has fewer points than ts 2 , the extra points in ts 2 (which are not coincident with any point in ts 1 ) will be excluded from the analysis. Pearson's correlation coefficient, ρ, is used to express the absolute relationship between any two parameters, without an attempt to predict the future. The values of ρ and d E averaged over all sites (ρ andd E , respectively) are also used to express the overall results in the study area. A crosscorrelation function (CCF), expressed as Eq. (9), is used to find the approximate time lag between the ATI derivations and the in-situ measurements.
ðx t −xÞðy tþk −ȳÞ½k ¼ 0; 1; : : : ; ðn − 1Þ; (9) where n is the total number of coincident observations in both series, t is time, k is lag, andx and y are the mean values of the x and y input series, respectively. 27

Short-Term Series: Apparent Thermal Inertia Results Versus Simulated and In-Situ Soil Moisture
Derived SM from the MODIS LST using the ATI method was compared with the in-situ and modeled SM at six test sites for a period of four months (August to November 2010). The accuracy of the WRF simulations is analyzed to find out if these simulations can be used to compensate for the gaps in ATI derivations. Comparison between the in-situ data collected from a singular point in space with the MODIS pixels (or the model's grid cells) was a major concern. The issue of pixel size versus point measurements, which is already dealt with in the literature, 28,29 is a well-known problem when it comes to defining the relationship between the in-situ and remotely sensed gridded datasets. This issue can be partly resolved through a careful selection of the measurement point closest to the pixel center over a large, flat, and homogeneous land area. Although all of our test sites were located on flat and relatively homogeneous landscape, the in-situ SM measurements only represented a fraction of the overlapping pixel or gridcell from the other two (i.e., ATI and WRF SM) datasets. As a result, a certain level of variation inside each pixel was inevitable. We tried to partly overcome this issue by choosing SM measurements at a 20-cm depth, which were considerably more homogeneous (see also Ref. 6).
As a result, a lower agreement between point measurements and gridded data is likely when a 1:1 correlation is applied. Therefore, more emphasis in the results will be given to the temporal profiles, and the correlation coefficient (ρ) values will be given only to provide a quantitative measure of the agreement between the two variables. ρ values were calculated between derived SM (using both ATI1 and ATI2), the WRF simulations, and the in-situ measurements ( Table 2). For most of the sites, ρ values between the WRF and the in-situ data have been positive and relatively higher than the ATI SM retrievals. Overall, ATI1 and ATI2 have correlated almost similarly with the in-situ measurements. Although these correlations showed almost similar results, temporal profiles of ATI1 and ATI2 SM retrievals needed to be compared to find out which approach presents a closer trend to the in-situ measurements and the WRF simulations. Therefore, temporal profiles of SM derived from both functions, simulated by the WRF model and measured on the ground, were overlaid on daily rainfall data (Fig. 3). The WRF simulations have shown a close match with the in-situ measurements, and rainfall effects are detected. The declining seasonal trend in SM amount, as observed on the in-situ profile, is also detected by the model but is less pronounced. Since ATI1 uses only two MODIS LST observations, there were more points possible to be calculated by this function. On the other hand, ATI2 needs four MODIS LST observations per day; thus, fewer points were possible to be modeled by this function. This is the reason that at some sites ATI2 has shown a weaker correlation (Table 2). However, it appears that the ATI2 profile is closer to the in-situ measurements than ATI1. Additionally, the declining seasonal trend is detected by the ATI2 function but not by ATI1.
There are more points in the later part of the analyzed period both in the ATI1 and ATI2 time series. This is due to the higher number of cloudy days in winter than in later spring. On the other hand, the model simulations have shown a larger offset with the in-situ measurements during the later part of the analyzed period. These results show that the ATI method is more useful in a less cloudy and warm season, unlike the WRF model, which performs better in a colder and wet season of the year. However, a longer period of analysis is required to confirm this finding, which is presented in the next section.

Long-Term Analysis: Apparent Thermal Inertia Results versus In-Situ Soil Moisture
The two ATI calculation functions, ATI1 and ATI2, were applied to derive the near-surface SM from the MODIS LST product for a period of three years (January 01, 2010, to December 31, 2012). Unlike the short-term analysis, this section will focus only on ATI SM derivations in comparison with the in-situ measurements and rainfall data.   Correlations among ATI1, ATI2, and the in-situ SM time series at six sites were calculated (Table 3). Because the in-situ measurements were made at a 20-cm depth below the surface, time lags were applied on these measurements before correlating with the derived SM time series. To discover the best time lag, the two time series (i.e, in-situ measurements versus ATI1 or ATI2 SM) were crosscorrelated with different time lags where every step in the lags was equal to one day. The best agreement on most sites was achieved when a time lag of one day was applied on the in-situ measurements. In other words, derived SM from the instantaneously observed LST data agreed best with the in-situ SM measured (at a 20-cm depth) one day later. The slope of change in correlations due to time lags, however, was very gradual with changes of less than 0.03 in ρ values. Thus, the actual time lag may have been slightly more or less than one day, which could be identified if data with a finer temporal resolution (such as hourly data) were available. However, retrieval of SM from the MODIS LST using the ATI method was only possible on a daily basis. Results showed that the correlations from both functions are relatively similar; however, ATI2 provided a better agreement with the in-situ measurements at more sites ( Table 3).
As discussed before, although these correlations offer a quantitative measure of how both functions relate to the in-situ measurements, they do not provide information about temporal trends of the derived SM with respect to the ground measurements. Therefore, temporal profiles of the derived SM time series at all sites were compared with the in-situ measurements and overlaid on rainfall data. The results from only two sites, Rangiora and Winchmore, are shown (Fig. 4). The most striking difference between ATI1 and ATI2 SM retrievals, as shown in the long-term series (Fig. 4), is that the latter has detected the annual SM trends better than ATI1. As found in the short-term analysis, both functions have shown a better performance in a dry season. The overall bias from the in-situ measurements was higher during winter compared to summer. To show this quantitatively, the SM time series from both functions were separated based on winter and summer seasons and then were subtracted by the corresponding in-situ measurements at each site. The mean error (ME) values at each site for both seasons were then calculated (Table 3). These values showed that in summer, both functions have a small bias from the in-situ measurements. The ME values from the in-situ measurements during the wet season were considerably higher, especially in ATI1 results (Table 3). Although rainfall effects are not clearly distinguished in this temporal scale, sudden spikes in ATI2 SM are suspected to be due to rainfall. On the other hand, the ATI1 series showed a relatively smooth temporal trend with only a few sudden spikes, also suspected to be due to rainfall. Another difference between ATI1 and ATI2 SM retrievals was the long-term offset between the two and the in-situ measurements. Although the ATI2 temporal trend overlaps the in-situ series in summer, both functions have shown a negative offset from the in-situ measurements for most of the three year period. To show this quantitatively, d E was applied. The time series of SM from ATI1 showed a larger difference with the in-situ time series than ATI2 at all sites (Table 3). It must be mentioned, however, that these d E values were calculated based on the absolute SM amount (m 3 · m −3 ) for the whole period; therefore, the one time series which had more points, tends to show more difference. As mentioned before, ATI2 had less points (due to the need for four daily LST observations) than ATI1 (which only needed two daily LST values). Therefore, the larger difference from ATI1 was partly due to more data points it its series. Since the ATI2 function uses both MODIS-Aqua and MODIS-Terra daily LST values, the DTA employed by this function contains more details. As a result, the temporal profile of SM derivations from this function is able to provide more details about SM variability in the long term. On the other hand, since four LST observations are required by ATI2, any missing LST has resulted in no SM retrieval for that day. This effect has been less severe on ATI1 as the possibility of having two cloudless LST per day is higher than four.
We also noticed that the prominent spikes due to rainfall events in the ATI series (especially ATI2) have occurred with a slight lag after the rainfall events. This was also shown in the correlations above (Table 3), where the best agreements were achieved when a time lag of one day was added to the in-situ SM measurements. It indicates that the changes in LST due to rainfall, as colder temperatures and smaller diurnal amplitudes, have been recorded by MODIS shortly after rainy days. Any changes in LST during the actual rain events (as a result of higher soil water content) are not available in the LST dataset due to cloud effects. As mentioned above, one way to fill these gaps in the ATI SM time series is to use the model simulations.

Discussion
Although satellite observations of SM based on TI are expected to capture the near-surface conditions, comparison with the in-situ SM measured at 2-5-cm depth turned out to be inconclusive in our investigation area. It is likely that these ground measurements have been affected by the local conditions and possibly direct solar radiation, hence, were highly variable within only a few meters' distance. This made a comparison not meaningful. Our analysis showed that SM measurements become relatively homogeneous over similar LC types in deeper soil (10 to 20-cm depth). Hence, to avoid the anomalies in the near-surface SM measurements, the authors used root-zone SM in the analysis after testing the time lag potentially required for moisture to reach a 20-cm depth in the soil.
The authors interpolated the 8-daily MODIS albedo product to daily for use in the ATI model. This approach may miss changes in surface albedo due to rainfall and cloudy conditions over the period of 8-day intervals. These changes, however, do not have a significant effect on LST-based SM derivations and can be cancelled out by the ATI model due to the following reasons. Rainfall affects the albedo of both wet and dry soil reducing the overall reflectance and surface temperature. Since we know that ATI is based on the difference between the maximum and minimum daily LSTs (Sec. 4), it can be assumed that rainfall will have a relatively similar effect on both temperatures. LST data are also not available when there is cloud cover, hence, no SM can be retrieved in cloudy conditions (see Fig. 3).
Our results showed that the agreement between the in-situ measurements at some sites, such as Leeston and Darfield, and the ATI SM retrievals were relatively poor in both short-term and long-term analyses. The Leeston site was close to Lake Ellesmere. Thus, it is suspected that the MODIS LST data over this site have been affected by water. This can be partly due to a spatial mismatch between the actual LST pixel overlying this site and the neighboring pixels over the lake. However, geometric mismatching in the MODIS LST grid had been already checked by overlaying this dataset on other spatial data from the study area (such as coastal boundaries and rivers) in a GIS (Geographic Information Systems) environment (Fig. 5). Proximity to mountains or local effects (such as irrigation) can be the reason for the poor correlations from the Darfield site. Irrigation is practiced widely during the summer to keep the grass growing in the farmlands across the study area. Irrigation effects may have caused higher SM in ATI estimations at those sites for the summer season, which can reduce the overall agreement with the ground measurements. Similarly, irrigation can cause anomalies in the in-situ data if the measurement site is not well located. It can cause significant divergence from the normal SM trend over a localized small area, which is not necessarily captured by the satellite observations or computed by the WRF model. Since the in-situ measurements were acquired from an already existing online database (Sec. 3), it was not possible for the authors to make sure that the measurements were not affected by irrigation. However, a quality check of the SM measurements based on rainfall data demonstrated that there were no external anomalies in the data except for the effects of rainfall. This can be checked in Figs. 3 and 4, where the only reasons for changes in SM appear to be due to rainfall and seasonal temperature variations.
The results also showed that the WRF model has agreed well with the in-situ measurements. Although this analysis was only conducted in the short-term period, these high correlations indicate that the WRF model can also be relied on for a longer period. As explained in Sec. 3.4, the initial conditions for the WRF coupled with the Noah LSM are based on the NCEP reanalysis data. The reanalysis data are produced operationally every 6 h by NCEP with a 40-km spatial resolution 12 based on the upper air and surface observational data collected from the measurements in the local weather stations via the global telecommunications system (GTS). The model also combines land surface parameters, such as the surface LC and vegetation fraction 12 to produce simulations. As a result, the ability of the model to incorporate observational data in the land surface parameterization schemes can be the reason for a good agreement with the ground measurements. However, there are two main drawbacks for a modeling system as opposed to a remotely sensed approach. First, the spatial resolution of the model is limited and the uncertainties increase when downscaling to a very fine resolution (≤1 km) is required. The second limitation is the computational time and cost for model simulations as well as the dependence of simulations on the availability of re-analysis data, which itself depends on the availability of measurements from permanent weather stations [Sec. 3(4)]. Satellite data, on the other hand, are available nearly real-time across the globe regardless of topography or distribution of weather stations and are readily available with a predefined spatial resolution. The WRF model coupled with the Noah LSM, therefore, can be a complementary solution to fill the gaps in SM retrievals from the MODIS LST data. An intercomparison of simulated and satellite-based SM by itself also can be useful for future research on both numerical modeling and/or remotely sensed retrieval of SM over moderate vegetation (mainly grass as in our study area) in New Zealand.
Although a simple 1:1 correlation helped us to find out which function offers a higher agreement with the in-situ SM measurements, comparison of ATI1 and ATI2 temporal profiles enabled us to figure out which function is better able to detect seasonal trends and rainfall effects. Temporal profiles also assisted us with the interpretation of the declining seasonal trends and rainfall effects in the WRF simulations. Overall, the ATI2 retrievals showed a better agreement with the in-situ measurements and rainfall events than the ATI1 outputs (Table 4). This indicates that using four daily LST observations from MODIS presents a better approximation of SM than the two (a minimum and a maximum) LST used by ATI1. This is because more often LST observations enable mapping the TI of the surface and interpretation of the day and night temperature difference. 30 However, the drawback of using four daily LST observations is that more uncertainty due to the satellite ground track variations is introduced into the ATI model. Additionally, using four LST observations by ATI2 results in more missing SM retrievals as opposed to ATI1, which needs only two daily LST values.
The lower correlations at some sites (Leeston and Darfield) and the missing SM values revealed that a remotely sensed approach for SM retrieval is limited to favorable weather conditions and suitable LC types. As the results showed, the accuracy of SM retrievals deteriorate over water, dense vegetation, and rugged terrain. Sensitivity of LST to SM differs for the canopy and the soil surface beneath the plants, and is much greater for bare soil than for canopies. 31 As a result, the accuracy of the ATI algorithm diminishes over dense vegetation. The retrievals are not possible under cloudy conditions. These limitations indicate that a remotely sensed method for SM retrieval can work only in regions with favorable conditions. Such a method, therefore, is not suitable for parts of New Zealand that have dense vegetation and rugged topography. Nevertheless, an intercomparison of simulated and retrieved SM over the moderately vegetated (mainly grass) LC types, as assessed in this study, can be useful for both modeling and satellite-based studies in the future.

Conclusion
Soil moisture derived from remotely sensed LST, simulated by a numerical model, and measured on the ground, was analyzed in this paper. The objective of the analysis was to understand the potential of the MODIS LST dataset for soil moisture retrieval using the ATI method and to compare the accuracy of two ATI functions based on the ground measurements. A land-atmospheric coupled model was also evaluated for potential gap filling of ATI soil moisture retrievals. Both ATI functions showed almost similar results in the short-term (four months) period, but the overall correlation between the ATI2 time series and the in-situ measurements was slightly higher (ρ ¼ 0.63) than ATI1 (ρ ¼ 0.62). At some sites (such as Leeston), both functions showed relatively poor correlations. It was discussed that the poor results at those sites were due to the effects of the nearby water bodies or mountains. The WRF simulations, on the other hand, showed relatively strong correlations at all sites (ρ ¼ 0.81). The model simulations also agreed well with the in-situ measurements in detection of rainfall effects and the general seasonal trend. Over the long-term analysis (three years), ATI2 showed slightly higher correlation (ρ ¼ 0.66) than ATI1 (ρ ¼ 0.63). Temporal profiles of the two functions showed a considerable offset from the in-situ time series in the long-term analysis; however, the overall bias in ATI2 retrievals was lower (d E ¼ 126) than ATI1 (d E ¼ 237), partly due to the fewer points in its series. To break down this overall bias into seasons, ME values during summer and winter for all sites were calculated. Both functions showed small biases from the in-situ measurements in summer (∼5 m 3 · m −3 volumetric SM), but considerably larger biases in winter with a slightly better result from ATI2 (ME ¼ −14 m 3 · m −3 ) compared to ATI1 (ME ¼ −14 m 3 · m −3 ). The ATI2 temporal profile was able to detect seasonal variations of SM better than ATI1. It was discussed that this is due to the more detailed DTA employed by ATI2 compared to the simple DTA of ATI1.
Results of this research indicate that the ATI2 function is more suitable for SM derivations in the study area. Since the MODIS LST product is available for more than 10 years and the mission is still being continued, long-term series of SM can be derived using the ATI2 function. Considering the good performance of the WRF model, the gaps due to cloud cover in ATI retrievals can be filled in by the model simulations.