The Importance of Subsurface Processes in Land Surface Modeling over a Temperate Region: An Analysis with SMAP, Cosmic Ray Neutron Sensing and Triple Collocation Analysis

: Land surface models (LSMs) simulate water and energy cycles at the atmosphere–soil interface, however, the physical processes in the subsurface are typically oversimpliﬁed and lateral water movement is neglected. Here, a cross-evaluation of land surface model results (with and without lateral ﬂow processes), the National Aeronautics and Space Administration (NASA) Soil Moisture Active/Passive (SMAP) mission soil moisture product, and cosmic-ray neutron sensor (CRNS) measurements is carried out over a temperate climate region with cropland and forests over western Germany. Besides a traditional land surface model (the Community Land Model (CLM) version 3.5), a coupled land surface-subsurface model (CLM-ParFlow) is applied. Compared to CLM stand-alone simulations, the coupled CLM-ParFlow model considered both vertical and lateral water movement. In addition to standard validation metrics, a triple collocation (TC) analysis has been performed to help understanding the random error variances of different soil moisture datasets. In this study, it is found that the three soil moisture datasets are consistent. The coupled and uncoupled model simulations were evaluated at CRNS sites and the coupled model simulations showed less bias than the CLM-standalone model ( − 0.02 cm 3 cm − 3 vs. 0.07 cm 3 cm − 3 ), similar random errors, but a slightly smaller correlation with the measurements (0.67 vs. 0.71). The TC-analysis showed that CLM-ParFlow reproduced better soil moisture dynamics than CLM stand alone and with a higher signal-to-noise ratio. This suggests that the representation of subsurface physics is of major importance in land surface modeling and that coupled land surface-subsurface modeling is of high interest.


Introduction
Soil moisture exerts an important control on the water and energy cycles in the atmosphere-land surface-subsurface continuum. Therefore, improving soil moisture estimation is beneficial for understanding the partitioning of water and energy fluxes. Soil moisture variability in the unsaturated zone is affected by water exchanges between the unsaturated zone, atmosphere, and groundwater. Studies have been focused on the effects of groundwater dynamics on land surface processes, showing the role of the groundwater in the water and energy cycles [1,2]. Results indicate that groundwater has little impact on soil moisture in deep groundwater regions, however, in districts with shallow groundwater-such as wetlands and river valleys-groundwater can become a major source of soil water [3][4][5]. The groundwater table depths and hydraulic gradients between saturated and unsaturated soils can cause capillary rise and make groundwater a constant water supply, leading to altered runoff and evaporation rates [6,7]. The land surface energy balance is affected by soil moisture states and land surface temperature [8][9][10].
Despite the importance of subsurface flow, most current land surface models (LSMs) neglect the lateral flow between grid cells or at sub-grid scales and only consider the water exchange in the vertical direction. Such models include the Variable Infiltration Capacity model (VIC), Interaction Soil Biosphere Atmosphere (ISBA) surface scheme, Noah Land Surface Model, and Community Land Surface Model (CLM). Recently, some works have considered the role of lateral subsurface flow and developed three dimensional hydrological models coupled with LSMs, such as CATHY (CATchment HYdrology), NoahMP [11], VIC-MD (MODFLOW) [12], and CLM-ParFlow [13]. These models aim at simulating the subsurface water and energy cycles more realistically than uncoupled models. However, model predictions are affected by errors given the uncertainty of the many required input parameters, e.g., atmospheric forcing, soil texture, and vegetation properties [14]. Precise soil hydraulic parameter data and land cover type information are hard to obtain as most areas lack sufficient land surveys. Also, the model structure and further assumptions influence the simulation performance. A large number of parameters used in land surface models are hard-coded as constants, although they are calculated by linear regression from preliminary studies using a limited amount of data and known to be uncertain [15].
Soil moisture can be measured via automated techniques such as gravimetric methods, nuclear techniques (such as neutron scattering, Gamma ray attenuation), electromagnetic methods (e.g., the time domain reflectometry (TDR) and the frequency domain reflectometry (FDR)), tensiometers and hygrometry [16][17][18]. However, most sensors monitor soil moisture at point scales (radius less than 1 m) only. As soil moisture is very heterogeneous in space and time, one usually needs to collect data from multiple locations in a specific area [19]. To reduce the scale gap to remote sensing products and modeling results and to obtain area-averaged soil moisture, a new technology has emerged with the Cosmic Ray Neutron Sensor (CRNS). It measures neutron count intensity and determines soil moisture in a non-invasive and continuous way [20]. The omnipresent cosmic radiation produces neutrons that interact with atmosphere and ground. These secondary neutrons include fast neutrons, that are generated by collisions between high-energy neutrons and nuclei. Fast neutrons are easily moderated by hydrogen atoms. As soil is the main source of hydrogen, thus, the variations of fast neutrons are strongly related to soil moisture changes. The process of moderation can be captured and counted by cosmic ray neutron sensors [21,22]. The large spatial footprint makes it suitable for agricultural water resources management and remote sensing product validation [23].
Two innovative satellite missions that include L-band (1200-1400 MHz) passive microwave systems have already been launched, including SMOS (Soil Moisture and Ocean Salinity mission, launched in November 2009) [24] and SMAP (Soil Moisture Active Passive mission, launched in January 2015 and starting operations in April 2015) [25]. Both SMOS and SMAP aim to measure soil moisture globally. After SMAP radar stopped operation in July 2015, an enhanced product was developed to give high-resolution observations at 9 km resolution. This enhanced product was evaluated by comparing it with long-term in situ soil moisture data [26,27]. It was found that the average ubRMSE (unbiased Root Mean Square Error) in L2_SM_P_E (Level 2 Enhanced Passive Soil Moisture) product and L3_SM_P_E (Level 3 Enhanced Passive Soil Moisture) product are between −0.040~0.055 cm 3 cm −3 [27][28][29], which barely meets the accuracy requirements of the SMAP mission. A number of previous studies have compared SMAP products and/or model simulations with in-situ observations. EI Hajj et al. [30] evaluated SMAP soil moisture products at sites in Southwestern France and found that the average bias over stations was about −0.032 cm 3 cm −3 , indicating SMAP moderately underestimates the soil moisture compared to in situ observations. Walker et al. [31] compared SMAP soil moisture with validation sites in the South Fork River watershed in Iowa, U.S., and found that the bias could be up to −0.04 cm 3 cm− 3 in early-spring and late fall and improve to −0.02 cm 3 cm −3 in the summer time. At global scale, a study indicated that SMAP shows a dry bias [32]. Recently, a new COSMOS network is developed to provide a unified, standardized, publicly available, traceable, and objective validation procedure that is operational in ISMN (https://www.geo.tuwien.ac.at/insitu/data_viewer/ accessed on 1 August 2021) and the QA4SM online validation service for soil moisture products (https://qa4sm.eu/ accessed on 1 August 2021) [33]. Here, soil moisture from ERA5, ERA5-Land, and GLDAS are provided for validation, but subsurface processes are not considered. As mentioned, SMAP was only compared to soil moisture simulated by stand-alone land surface models that have a rather simple subsurface structure and do not consider the role of groundwater and lateral flows [34][35][36], which could lead to systematic deviations. Studies show that soil moisture is overestimated when the model neglects the impact of topography [37,38]. Also, previous studies focused mostly on large scales and simulations at a coarser spatial resolution [23,39,40].
To find relative error estimates for different products, the Triple Collocation (TC) method is used, which is an error magnitude estimation approach for intercomparison among three or more independent observed or modeled datasets [41]. The method was firstly used for ocean wind studies and then developed and widely used in other areas, such as land surface hydrology [42][43][44][45][46][47]. It has been proven to be a useful tool to understand the random error variances of remote sensing time series. It assumes one of these datasets as reference to relative rescaling, and further assumes truth-error orthogonality and zero error cross correlation between datasets to obtain a bias-free TC analysis [46][47][48]. Here, we present simulated soil moisture by CLM and the CLM-ParFlow coupled model and compared it with CRNS observations and SMAP enhanced soil moisture datasets using TA method. Compared to the CLM model, CLM-ParFlow considers three-dimensional water flow in the subsurface (soil and aquifer) and a two-dimensional overland flow module. It is investigated whether a better subsurface representation can improve soil moisture estimates. The research area has various land use types and various pedological, hydrological, and hydrogeological site conditions were observed. These small-scale hydrological processes can provide insights for large scale modeling. This research aims to understand the performance and limitations of the subsurface role in land surface model simulations, which might recall the need to consider hydrogeology (including complex 3D geology and critical parameters like hydraulic conductivity and storage coefficients) in soil-vegetationatmosphere processes, and to obtain finally more accurate soil moisture and groundwater level data.

Study Area
The study area is located in central Europe encompassing parts of North Rhine-Westphalia and Rhineland-Palatinate in western Germany and parts of Belgium, the Netherlands and Luxemburg, covering an area of 150 × 150 km ( Figure 1). This region has a sub-Atlantic oceanic climate. Summers are mild, while winters are humid and relatively mild. The average monthly temperatures are highest in July (18 • C) and lowest in January (3 • C). There is a large spatial variability in precipitation due to topography. In the rather flat Lower Rhine area in the North of the study region, the yearly precipitation is between 600-900 mm [49]. In the southern low mountain ranges the average yearly precipitation is locally 1600 mm in the Bergisches Land (South East) and 1300 mm in the Eifel (South West) [50]. The rainfall is frequent and evenly distributed over the seasons. The elevation in this area ranges from the plains at around 14 m a.s.l. to the mountainous areas at around 735 m. The land use is a mixture of agriculture area, forests, urban and rural areas, wat grassland, industry, and mining [2]. The dominant land use is agriculture, covering m than 60% of the area (shown in Figure 2). The major crops are winter wheat, maize, a sugar beet. Forests cover nearly 20 % of this area and are mainly located in the south, i Eifel, Bergisches Land, and Sauerland. The dominant soil textures are sandy loam, loa and clay loam and our simulations are based on the FAO/UNESCO Soil Map [51]. San soils with low water holding capacities are mainly located in the Northwest region. This research area has three solid rock areas of regional importance [52]. The larg part is Palaeozoic rocks located in the southern part and east of the river Rhine, hav the dominant aquifer typology 'schist and shales' that includes folded and partly me morphosed clastic sedimentary rocks. This part has a low hydraulic conductivity. T north-eastern part is mainly occupied by unconsolidated rocks of Cenozoic era, includ alternating sequences of clastic sedimentary rocks horizontally. The hydraulic conduct ity varies a lot depending on the combination system. In the western part, the domin aquifer rocks are mostly consolidated, weakly permeable material from Mesozoic. The land use is a mixture of agriculture area, forests, urban and rural areas, water, grassland, industry, and mining [2]. The dominant land use is agriculture, covering more than 60% of the area (shown in Figure 2). The major crops are winter wheat, maize, and sugar beet. Forests cover nearly 20 % of this area and are mainly located in the south, i.e., Eifel, Bergisches Land, and Sauerland. The dominant soil textures are sandy loam, loam, and clay loam and our simulations are based on the FAO/UNESCO Soil Map [51]. Sandy soils with low water holding capacities are mainly located in the Northwest region. The land use is a mixture of agriculture area, forests, urban and rural areas, water, grassland, industry, and mining [2]. The dominant land use is agriculture, covering more than 60% of the area (shown in Figure 2). The major crops are winter wheat, maize, and sugar beet. Forests cover nearly 20 % of this area and are mainly located in the south, i.e., Eifel, Bergisches Land, and Sauerland. The dominant soil textures are sandy loam, loam, and clay loam and our simulations are based on the FAO/UNESCO Soil Map [51]. Sandy soils with low water holding capacities are mainly located in the Northwest region. This research area has three solid rock areas of regional importance [52]. The largest part is Palaeozoic rocks located in the southern part and east of the river Rhine, having the dominant aquifer typology 'schist and shales' that includes folded and partly metamorphosed clastic sedimentary rocks. This part has a low hydraulic conductivity. The north-eastern part is mainly occupied by unconsolidated rocks of Cenozoic era, including alternating sequences of clastic sedimentary rocks horizontally. The hydraulic conductivity varies a lot depending on the combination system. In the western part, the dominant aquifer rocks are mostly consolidated, weakly permeable material from Mesozoic.
Soil moisture validation activities have been performed in this area before based on in situ and CRNS data [23,[53][54][55][56] as well as simulation and data assimilation experiments [57][58][59][60]  This research area has three solid rock areas of regional importance [52]. The largest part is Palaeozoic rocks located in the southern part and east of the river Rhine, having the dominant aquifer typology 'schist and shales' that includes folded and partly metamorphosed clastic sedimentary rocks. This part has a low hydraulic conductivity. The north-eastern part is mainly occupied by unconsolidated rocks of Cenozoic era, including alternating sequences of clastic sedimentary rocks horizontally. The hydraulic conductivity varies a lot depending on the combination system. In the western part, the dominant aquifer rocks are mostly consolidated, weakly permeable material from Mesozoic.

Land Surface Modeling
Soil moisture was simulated by the land surface model CLM (community land model) and the coupled land surface-subsurface model CLM-ParFlow. The CLM model was developed by the National Center for Atmosphere Research (NCAR). In the CLM 3.5 model, the soil is discretized into 10 unevenly distributed soil layers (see Table 1). Soil hydraulic properties are estimated internally from soil texture (sand fraction and clay fraction) using pedo-transfer functions according to Clapp and Hornberger [61] and Cosby et al. [62]. A simplified Richards equation is used in CLM to calculate the vertical water movement in the unsaturated zone ∂θ ∂t where ∂θ ∂t is the soil moisture (cm 3 cm −3 ) change over time, K represents the saturated hydraulic conductivity (m/s), ϕ is the pressure head (unit is defined as length (L)), z as the vertical coordinate (L) and Q is the source/sink term (i.e., the soil water removed due to evaporation). A main limitation of CLM 3.5 is that lateral flows are not considered, and groundwater is not represented, although groundwater can strongly influence soil moisture conditions. ParFlow [63] is a numerical, integrated hydrological model that simulates subsurface groundwater flow and water flow in soils, as well as overland flow. Both retention and relative permeability curves are represented by the van Genuchten relationships [64]. ParFlow does not include land surface process (e.g., evaporation), nor does it have a parameterization scheme for frozen soil and ice processes. In addition to the DEM dataset used in CLM, the topographic slopes need to be specified for ParFlow. When coupled with CLM, ParFlow replaces the one-dimensional CLM soil moisture characterization by the three-dimensional approach in CLM-ParFlow, considering the redistribution of soil moisture, and integrating vertical and lateral flow of groundwater and surface water.
In this equation, S s is the specific storage coefficient (L −1 ), S w is the relative saturation, φ is the porosity, k r is the relative permeability. The subsurface is discretized into 30 layers, with 10 vertically layers near the surface (2-100 cm) and 20 constant levels (135 cm depth) that reach up to 30 m below the surface. The physically based coupled model (CLM-ParFlow) can better simulate the role of groundwater in terrestrial systems, and the interaction between surface water and subsurface [13,63].
To represent the high spatial heterogeneity of the land surface, the simulation domain was discretized into grid cells of 500 × 500 m. The plant functional types (PFTs) were based on MODIS land cover data. The hilly areas are mostly covered by broadleaf forest and needleleaf forest, while the other flat regions are covered by crops and urban areas. Soil texture information was taken from the FAO/UNESCO Soil Map [51] with the scale of 1: 5,000,000 (see Figure 2). Most of the model domain is dominated by clay-loam (35% clay, 35% sand, 30% silt). Sandy loam (10% clay, 65% sand, 25% silt) and loam (20% clay, 40% sand, 40% silt) are dominant in the north-western part. Clay (45% clay, 15% sand, 40% silt) is found in the northwestern corner of the domain. The van Genuchten water retention Remote Sens. 2021, 13, 3068 6 of 22 parameters and the hydraulic conductivity used in CLM-ParFlow are calculated by the Rosetta pedo-transfer function [65], summarized in Table 2. To drive the model, the high-resolution reanalysis dataset (COSMO-REA6) [66] was used as meteorological forcing. This dataset covers the period 1995-2020 and is continuously supported by DWD (Deutscher Wetterdienst, German Meteorological Service). It uses ERA-Interim data as a boundary condition and is generated by assimilating observed meteorological data into the atmospheric model COSMO (Consortium for Small-scale Modeling) [67]. The dataset comprises air temperature, precipitation, humidity, incoming shortwave and longwave radiation. A two-year spin-up period was applied for CLM and the initial conditions for CLM-ParFlow were taken from previous studies [2,59]. Both simulations with CLM and CLM-ParFlow started from a near equilibrium condition. Thus, different spin up treatments do not have influence on results. In total, a period of two years (2017-2018) was simulated with a time step of 3600 s. In case of convergence issues, the time steps are reduced until convergence can be achieved.
Simulated soil moisture for the upper 5 cm layer (SM 5cm ) and the upper 20 cm (SM 20cm ) were estimated by linearly combining simulated output for different model layers (H20SOI i , where i is the index denoting the soil layer).

CRNS Observations
The CRNS is an emerging technology to monitor soil moisture at the intermediate scale [20,22]. The measured neutron count intensity provides an estimate of soil moisture content for a radius of around 240 m, at sea level and dry bare soil conditions. The radius is a function of air density, air humidity and vegetation density [68]. The penetration depth of the CRNS measurements varies from 15 cm (wet soils) to 55 cm (dry soils) [69]. The neutron count intensity is mainly sensitive to the number of hydrogen atoms in the soil, but is also influenced by changes in atmospheric pressure, vapor pressure, and incoming cosmic radiation. These factors are considered in the standard correction process [69,70].
Several studies have been conducted to investigate the accuracy of the CRNS measurements and found that CRNS provides reliable soil moisture estimates when calibrated properly [71][72][73]. Bogena [74] found that even for a densely vegetated and wet site, the RMSE of daily soil moisture estimated by CRNS is only 0.03 cm 3 cm −3 . In our work, 13 CRNS stations are used to evaluate SMAP and soil moisture model products (see Figure 1). The datasets are collected in the context of the TERENO project [75] and passed quality assurance procedures. We acknowledge that the effective depth of CRNS is dependent on soil moisture, and also on the depth of the calibration dataset. Our CRNS calibration and validation used individual support volumes of samples from 5, 20, and 50 cm based on the gravimetric method. We calculate the penetration depth based on previous study [76] and found that the effective penetration depth is mostly between 15-30 cm in our research area and period. The mean and median of penetration depth are 20.3 cm and 18.6 cm respectively. Hence, we assumed that the CRNS has an effective penetration depth of 20 cm. Table 3 provides more information about the CRNS stations. Table 3. Coordinates, altitude (m), average annual precipitation (mm y −1 ), land use type information [60] for 13 sites, adapted with permission from ref. [60]. Copyright 2017 by Roland Baa.tz.

SMAP Enhanced Soil Moisture Product
SMAP provides soil moisture observations of the top 5 cm of the soil and thaw/freeze states derived from the passive microwave brightness temperature (BT). BT is recorded by a conically-scanning antenna beam at L-band with a 40 • incidence angle. This results in a −3 dB antenna footprint of 40 km. To enhance the resolution of the typically 36 km SMAP radiometer data posting, the Backus-Gilbert optimal interpolation technique is used to interpolate the multiple scans of a single location. It makes most use of the available information and provides a better representation of the original data [26].
In this study the L3_SM_P_E product (version 4) was used, which provides soil moisture on a 9 km EASE2 (updated Equal-Area Scalable Earth-2) grid (National Snow and Ice Data Center NSIDC, https://nsidc.org/data/smap/smap-data.html, accessed on 1 August 2021). The soil moisture baseline retrieval algorithm in L3_SM_P_E product is performed by the vertical polarization single channel algorithm (SCA-V) (https://nsidc. org/data/smap/technical-references, accessed on 1 August 2021). The L3_SM_P_E product is provided in the form of global daily datasets, including soil moisture measured for the 6:00 a.m. (descending) and 6:00 p.m. (ascending) orbit. Here, the soil moisture daily value is calculated by taking the average of the two datasets. To eliminate the non-high-quality pixels, the surface and quality flags are used (retrieval_qual_flag and surface_flag).

Data Processing
Both 2017 (normal year) and 2018 (dry year) were selected to be evaluated for all datasets. In order to avoid unreliable soil moisture observations during frozen conditions and snow cover, the winter period (December, January, and February) was excluded.
The model simulations and SMAP product are sampled onto the SMAP grid by nearest neighbor (NN) search. Compared to area-wide spaceborne observation and model simulation results, the CRNS stations are quite sparse. As the spatial coverage of measurements by a CRNS is close to the model grid size in this work, CRNS observations are compared to a complete cell of the CLM, CLM-ParFlow, and SMAP grid containing the coordinates of the CRNS stations. The SMAP product and the model grids are at different spatial resolutions (9 km vs. 500 m) and comparison between SMAP soil moisture and modeled soil moisture is made at both resolutions. The fine resolution (model grid) provides a detailed assessment of the spatial variability, and the coarse resolution (SMAP grid) gives a smoothed representation excluding local noise. The comparison at the 9 km resolution is made by taking the arithmetic average of simulated soil water content data on the 500 m grid to get modeled values at the 9 km resolution (upscaling). The comparison at the 500 m resolution is made by downscaling the SMAP data, using the nearest source to destination to remap from 9 km resolution to 0.5 km resolution. In this case, for each model grid cell, it takes the value from the nearest SMAP grid cell. This is done by ESMF (Earth System Modeling Framework) regridding function in NCL (NCAR Command Language).

Standard Evaluation Metrics
Statistical performance was evaluated according to Good Practices Guidelines [19,77], including bias, Root Mean Square Difference (RMSD), and unbiased Root Mean Square Error (ubRMSD). The bias is given by where X i represents a simulated or remotely sensed product, and X i,re f is the referenced soil moisture dataset. The sample size is n. The Root Mean Square Deviation (RMSD) is given by and the unbiased root mean squared difference (ubRMSD) is The correlation between different soil moisture datasets (X and Y) was calculated using the Pearson correlation coefficient (r X,Y ). Here, the focus is on the dynamics of the different time series rather than absolute soil moisture values, considering the fundamental scale mismatch between the CRNS, SMAP, and LSM model grid, which cannot be adequately considered without complex scaling functions. The ubRMSD and r X,Y screen out the influence of bias between different amplitudes of soil moisture variation [78].
where cov(X, Y) represents the covariance between two datasets, σ X and σ Y are the standard deviations for dataset X and Y.

Triple Collocation
The method is based on the assumption that a measurement system can be understood as being composed of an additive systematic error (α i ), multiplicative systematic error (β i ), additive zero-mean random error (ε i ), and θ (the underlying truth value) [45] X, Y, and Z denote the time series of soil moisture products to be compared. Here, the assumption is that the errors are uncorrelated with each other (cov ε i , ε j = 0, i = j) and with θ (cov(ε i , θ) = 0). It applies for both reference datasets and soil moisture products to be evaluated. Given that all observation data has errors, in this study, we assumed that X is the reliable CRNS dataset which is well calibrated and without multiplicative (second order) error, Z is the SMAP L3 enhanced product, and Y denotes the soil moisture in the upper 5 cm of CLM and CLM-ParFlow simulations. The β x is set to 1, and other scaling factors are calculated according equations 10-12 to eliminate the scale differences of different products. This assumption was also made in other studies [45,77,79].
The scaling factor can be regarded as a form of regression, by taking a third variable as a tool to resolve the relationship between the two variables [80]. It can be used to describe the soil moisture sensitivities. The unscaled error variances which represent the variance of the scaled white noise of each product can be solved by In addition, the signal-to-noise ratio (SNRs) is calculated, which provides a more clear representation of the ratio between soil moisture and uncertainty magnitude [77]. It is a combination of several information sources, including the sensitivity of the measurement system, the variability of the true value (θ 2 ) and the variability of the random error (ε 2 ) [44]. SNRs is not expressed as normalized between 0 and 1, but is often calculated and linearized by converting into decibel (dB) units Table 4 shows that the SMAP L3_SM_P_E product and CRNS observations have in general a relatively good agreement. The detailed time series of soil moisture from SMAP, CRNS and LSMs at 13 sites are provided in Appendix A ( Figure A1). The SMAP product and CRNS show on average not large systematic differences, while previous studies [23,39] found that SMAP tends to underestimate soil moisture compared to CRNS measurements at most sites. In terms of average r value, the SMAP is relatively well correlated to CRNS, ranging from 0.653 to 0.825, except for Ruraue station (r = 0.452). The Ruraue station is located near the Rur river, and part of the closest SMAP pixels contain some amount of open fresh water. The presence of open water introduces a soil moisture bias due to the lower brightness temperature for the grid cell. This may partly explain the low r for the Ruraue station. However, it should be noted that local differences are large. One reason is that the CRNS footprint is still much smaller than a 9 km SMAP-pixel. This spatial mismatch could lead to differing soil hydraulic parameters and land cover between CRNS and SMAP observations (see Tables 3 and 5) causing lower correlation between the two soil moisture datasets. It is also important to consider that the satellite observations are negatively impacted by high vegetation density, topography, frozen soil, snow cover, and volume scattering effect in case of low soil moisture content. The retrieval under dense forest is challenging or impossible since the recorded signal (brightness temperature) originates to a large extend from canopy instead of the soil microwave emissions [81,82]. When vegetation density increases, the impact of soil moisture on changes in ground emissivity becomes invisible, hence, the contribution of the ground is less than from the canopy. A large bias is observed in Wildenrath, where the land cover type is needleleaf (forest). Moreover, topography as well as surface roughness increases the surface area and alters the total microwave emission. In addition, it is also found that areas with complex topography are prone to shadowing and adjacency effects [83,84].

Agreement between Spaceborne and In Situ Observations
It should be noted that the soil organic matter content is high at the Wüstebach site, and there is a large difference of the bulk density between ancillary data (1.3 g cm −3 ) and actual observed value (0.83 g cm −3 ). In the single channel algorithm used by SMAP to retrieve soil moisture, the dielectric mixing model plays an important role in describing the relationship between soil moisture and microwave emissivity. A recent study has found that high levels of organic matter decrease the microwave effective dielectric constant and therefore cause higher brightness temperature for a particular soil moisture content [85].

Comparison of Model Simulation and CRNS Measurements
The 500 m modeled soil moisture was compared with CRNS measurements, assuming a measurement depth of 20 cm. Bias, RMSD, ubRMSD, and correlation coefficient were calculated for 13 stations and over three seasons (see Table 6 and Figure 3). The CLM simulations tend to overestimate soil moisture with a bias of 0.070 cm 3 cm −3 and CLM-ParFlow has a slight dry bias of −0.021 cm 3 cm −3 . The large wet bias of the models for Wildenrath is most likely related to soil texture. Soil moisture decreases faster in sand than in finer textured soil. The sand content is high (up to 65%), and the models seem to have a too low hydraulic conductivity for this site. CLM-ParFlow shows a large bias at sites located at a higher elevation, such as Wüstebach, Rollesbroich, and Schoeneseiffen sites.   Figure 4 shows daily mean soil moisture content in the upper soil layer (5 cm) over the research area for CLM and CLM-ParFlow simulations compared with the SMAP L3_SM_P_E product, together with daily precipitation time series. The SMAP product and model simulations have different soil wetting and drying dynamics after rainfall. The near-surface soil is sensitive to precipitation because of intensive positive vertical water gradients, but tends to dry quickly after water infiltration related to evaporation [86]. Although the representation of L-band sensing depth at 5 cm has been using for remote sensing validations [26,87,88]; however, the L-band sensing depths are affected by soil moisture, and the penetration depth can be shallower when soil moisture is high [89]. The penetration depths and vertical soil moisture gradients lead to different drying behavior. The simulations match soil moisture from SMAP well during most of the year. Generally, CLM overestimates the soil water content during summer. Compared to SMAP, CLM-ParFlow simulations show a very small wet bias of only 0.004 cm 3 /cm 3 while CLM has a larger wet bias of 0.065 cm 3 cm −3 . The RMSD for CLM and CLM-ParFlow simulations are 0.085 cm 3 cm −3 and 0.045 cm 3 cm −3 respectively.  After overlooking systematic bias, ubRMSD for the CLM stand-alone model and the coupled model are not significantly different and both below 0.06 cm 3 cm −3 , which indicates that both simulations could produce reasonable results. The correlation between CRNS observations and land surface simulations ranges from 0.576 to 0.821 for CLM stand-alone and 0.365 to 0.805 for CLM-ParFlow. Both CLM and CLM-ParFlow simulations show a strong correlation with in situ data, suggesting that the soil moisture dynamics at 500 m scale can be relatively well captured by both models. Figure 4 shows daily mean soil moisture content in the upper soil layer (5 cm) over the research area for CLM and CLM-ParFlow simulations compared with the SMAP L3_SM_P_E product, together with daily precipitation time series. The SMAP product and model simulations have different soil wetting and drying dynamics after rainfall. The near-surface soil is sensitive to precipitation because of intensive positive vertical water gradients, but tends to dry quickly after water infiltration related to evaporation [86]. Although the representation of L-band sensing depth at 5 cm has been using for remote sensing validations [26,87,88]; however, the L-band sensing depths are affected by soil moisture, and the penetration depth can be shallower when soil moisture is high [89]. The penetration depths and vertical soil moisture gradients lead to different drying behavior. The simulations match soil moisture from SMAP well during most of the year. Generally, CLM overestimates the soil water content during summer. Compared to SMAP, CLM-ParFlow simulations show a very small wet bias of only 0.004 cm 3 /cm 3 while CLM has a larger wet bias of 0.065 cm 3 cm −3 . The RMSD for CLM and CLM-ParFlow simulations are 0.085 cm 3 cm −3 and 0.045 cm 3 cm −3 respectively. ture, and the penetration depth can be shallower when soil moisture is high [89]. The penetration depths and vertical soil moisture gradients lead to different drying behavior. The simulations match soil moisture from SMAP well during most of the year. Generally, CLM overestimates the soil water content during summer. Compared to SMAP, CLM-ParFlow simulations show a very small wet bias of only 0.004 cm 3 /cm 3 while CLM has a larger wet bias of 0.065 cm 3 cm −3 . The RMSD for CLM and CLM-ParFlow simulations are 0.085 cm 3 cm −3 and 0.045 cm 3 cm −3 respectively.  Spatial maps of performance indices are given in Figure 5 and illustrate the differences between model predictions and observations. For the maps at 500 m resolution, as the resolution of the simulations is finer than the satellite measurement, large differences between simulation and measurements occur in valleys and rivers regions because these areas are not well covered by the coarser resolution of the SMAP satellite.

Temporal and Spatial Correlation between Model Simulations and the SMAP L3_SM_P_E Product
Remote Sens. 2021, 13, 3068 13 of Spatial maps of performance indices are given in Figure 5 and illustrate the diff ences between model predictions and observations. For the maps at 500 m resolution, the resolution of the simulations is finer than the satellite measurement, large differenc between simulation and measurements occur in valleys and rivers regions because the areas are not well covered by the coarser resolution of the SMAP satellite. In most of the area, CLM has a higher soil water content than SMAP. The CLM-P Flow model captures the spatial variability of soil water content and shows the influen of the river network. The soil moisture is close to saturation in the river valley. Meanwh the soil is drier in hilly areas. This is related to the difference in subsurface process rep sentation between CLM and CLM-ParFlow. In Parflow, the Richards' equation is used calculate 3D subsurface water flow, including both vertical and lateral water moveme which includes the lateral groundwater flow which moves water from the hilly areas wards the river valleys. In addition, lateral flow by streams and rivers is also modele The water flow convergence process (i.e., the lateral redistribution of water via strea and aquifers from hills to the river valleys and lowlands) can be better captured by CL ParFlow than by CLM, which just considers vertical water flow. In the northern flat a valley areas, the differences of soil moisture between CLM and CLM-ParFlow is small where precipitation and infiltration excess is low, thus the later water flow redistributi processes have smaller impact.
Further minor discrepancies between CLM and CLM-ParFlow simulations are lated to the different estimation of soil hydraulic parameters (different pedo-transfer fun tions). Several areas where CLM shows larger deviations with respect to SMAP soil mo ture coincide with loam regions of the soil texture map, which indicates that for loa texture CLM is probably too wet, which could be related to the default pedo-transfer fun In most of the area, CLM has a higher soil water content than SMAP. The CLM-ParFlow model captures the spatial variability of soil water content and shows the influence of the river network. The soil moisture is close to saturation in the river valley. Meanwhile, the soil is drier in hilly areas. This is related to the difference in subsurface process representation between CLM and CLM-ParFlow. In Parflow, the Richards' equation is used to calculate 3D subsurface water flow, including both vertical and lateral water movement which includes the lateral groundwater flow which moves water from the hilly areas towards the river valleys. In addition, lateral flow by streams and rivers is also modeled. The water flow convergence process (i.e., the lateral redistribution of water via streams and aquifers from hills to the river valleys and lowlands) can be better captured by CLM-ParFlow than by CLM, which just considers vertical water flow. In the northern flat and valley areas, the differences of soil moisture between CLM and CLM-ParFlow is smaller, where precipitation and infiltration excess is low, thus the later water flow redistribution processes have smaller impact.
Further minor discrepancies between CLM and CLM-ParFlow simulations are related to the different estimation of soil hydraulic parameters (different pedo-transfer functions). Several areas where CLM shows larger deviations with respect to SMAP soil moisture coincide with loam regions of the soil texture map, which indicates that for loam texture CLM is probably too wet, which could be related to the default pedo-transfer function in CLM, which might underestimate hydraulic conductivity for loamy soil types. The spatial distribution of both bias and ubRMSD also show clear differences between the northern flat area (with larger bias) and the southern hilly area (smaller bias). This is probably related to the fact that precipitation is high in the hills and both CLM and SMAP have soil moisture values close to field capacity. Concerning CLM-ParFlow, the differences with SMAP are larger in the Rhine valley, ParFlow overestimates the influence of the river on soil moisture values in the areas next to the stream. Notice that the CLM-ParFlow model is not calibrated, and that river-groundwater interaction could be closer to real conditions by adjusting for example riverbed hydraulic conductivities. Both bias and ubRMSD values show that the discrepancies between model simulations and SMAP data are smallest in autumn for CLM simulations and CLM-ParFlow.

Triple Collocation
In general, the statistics of a time series triple is inherently unique, so that a comparison of different triples in a collocation analysis is not directly possible. However, the same SMAP product and CRNS data are used in this study to be compared to two different simulation results. That provides a large ability for a direct comparison. Table 7 shows the TC results for model simulations (as Y) and SMAP L3_SM_E_P product (as Z) compared to reference CRNS datasets (as X). Per definition, the scaling factor β X of the reference dataset is 1, while the SMAP product and model simulations are scaled to this reference product. β y and β z values larger than 1 indicate that the dynamic range of the datasets from the SMAP product and the model simulations is lower than that of CNRS soil moisture time series, and vice versa. For CLM stand-alone simulation, the average β y is 4.538 and β z is 1.370. For CLM-ParFlow, β y and β z are quite close (1.453 and 1.242 separately). The scaling factors larger than 1 indicate that both SMAP L3_SM_P_E and the land surface model simulations have underestimated the soil moisture dynamics at the CRNS-sites, and probably in the complete research area. Compared to CLM stand-alone simulations, the β y of the coupled model is closer to 1, indicating less need to scale. In terms of β z , it should be noted that the scaling factors for Wüstebach and Wildenrath are lower than for other stations, indicating that the dynamics range of retrieved soil moisture needs to be reduced to match the scale of CRNS observations and model simulation time series. For these two stations, some additional problems (see Section 3.1) seem to influence the soil moisture retrieval results. Although the soil moisture datasets used in the TC correspond to different depths (SMAP for the upper 5 cm and the other two datasets for the upper 20 cm), no obvious relationship between the penetration depth differences and soil moisture dynamics is detected. Table 7 also provides the absolute TC error standard deviation. The ranges for σ Y are 0.010 to 0.026 (CLM) and 0.028 to 0.080 (CLM-ParFlow) respectively. In general, the SMAP soil moisture product provides similar error standard deviations with σ Z ≤ 0.060 in both triples.  The linearized SNR value gives the ratio relationship between soil moisture and uncertainty magnitude. On average, the SNR Y of CLM and CLM-ParFlow are 0.145 and 1.740 dB respectively. The negative SNR Y values demonstrate that the random noise is larger than the soil moisture signal. In the CLM-ParFlow simulation, most sites show a better performance than the CLM model. As described before, Ruraue is along the river, where the soil moisture is sensitive to river parameterization in ParFlow. Heinsberg, Kall, and Kleinhau have large negative values, indicating that both models have large absolute errors. Overall, the CLM-ParFlow could provide more valuable results than the stand-alone model.

Effect of Lateral Water Flow on Soil Moisture
Previous research indicated that the lateral flow is important in land surface modeling, especially when the resolution is fine [90]. The spatial patterns of modeled soil moisture show that the CLM-ParFlow has wet grid cells at foothills and valleys. The soil moisture gradient is larger in the wet grid cells and surrounding drier grid cells. By taking account of lateral flow, soil moisture decreases in these wet areas due to lateral diffusion. Also, the lateral drainage driven by topographic gradient results in soil moisture redistribution. In view of previous studies [91], lateral flow is expected for steep hillocks, even if slight difference in soil texture between adjoining grid cells. Also, the accumulated runoff in ParFlow, generated by infiltration excess or saturation excess, can route or reinfiltrate, while some other traditional LSMs can remove excess water from modeled water cycle [92], CLM-ParFlow maintains high soil moisture in convergence areas. This confirms the importance of lateral subsurface flow on the hydrological cycle, especially in mountainous areas.

Conclusions
This study compared soil moisture data from cosmic ray neutron sensors (CRNS), passive microwave remote sensing (SMAP L3_SM_E_P product) and land surface model simulations by the Community Land Model (CLM, version 3.5) and the coupled land surface-subsurface model CLM-ParFlow over a 150 × 150 km region in western Germany. CLM-ParFlow can better capture the impact of groundwater on soil moisture than CLM as it has a more advanced subsurface physical process scheme. With this approach an analysis of the impact of the representation of subsurface processes in hydrological simulations of soil moisture was performed. The evaluation results can be summarized as follows: Over 13 CRNS sites, the SMAP L3_SM_E_P product shows a small bias of 0.005 cm 3 cm− 3 only compared to the CRNS observations. Nevertheless, local differences can be large (up to −0.120 cm 3 cm− 3 for the densely forested Wüstebach site) due to differing spatial resolution of the soil moisture products and errors in soil texture and land cover used to derive soil moisture from brightness temperature in the SMAP L3_SM_E_P product. Besides, the disturbing role of dense vegetation and complex topographic features influence the accuracy of the SMAP product. Overall, the unbiased root mean square error (ubRMSE) is around 0.056 cm 3 cm− 3 , indicating that SMAP L3_SM_E_P product could barely meet its mission requirement for this very heterogeneous and hilly region.
The comparison between CRNS and land surface simulations show that CLM has a wet bias (0.070 cm 3 cm −3 ) and CLM-ParFlow has a dry bias (−0.021 cm 3 cm −3 ). Local biases can be large, which might be related to the uncertainty in soil texture and hydraulic conductivity, inadequate pedotransfer functions and lack of consideration for soil bulk density in CLM model. In terms of ubRMSE, both CLM and CLM-ParFlow are below 0.06 cm 3 /cm 3 and compare well to CRNS observation dynamics. The SMAP product and CLM-ParFlow do not show a systematic difference in soil moisture, which is in contrast to most land surface models which are wetter than SMAP.
The triple collocation (TC) comparison implies that both CLM and CLM-ParFlow show similar noise levels with σ Z below 0.058. The scaling factor of CLM-ParFlow is less than a third of CLM stand-alone, indicating that the coupled model could perform better with respect to CRNS measurements. This is an important aspect for future data assimilation studies, as the typical adaptation of the soil moisture climatology of model and observation becomes less mandatory. The higher SNR (signal-to-noise ratio) value for the coupled model CLM-ParFlow also indicates it can provide more valuable results than the CLM stand-alone model.
It should be noted that the direct metrics (e.g., RMSE and r) do not show a clear better performance of the CLM-ParFlow model compared to the CLM stand-alone model. The TC method shows that the simulation has been improved when lateral subsurface dynamics is involved. Unlike typical performance metrics, where the assumption is that the reference data set is free of (random) errors, TC methods account for sensor and representativeness errors and can be considered more robust than conventional metrics and close to reality [45]. With conventional evaluation metrics, we focus on the dynamics of the different time series instead of the absolute soil moisture values, because there can be a systematic bias between CRNS and SMAP measurements, as well as model simulations, which is related to different underlying assumptions for the different measurement and simulation methods. This method is also used because it is standard in the land surface modeling literature and allows therefore an easier comparison with other papers [36,40,93,94].
In summary, the model structure is important for soil moisture modeling. Compared to CLM-ParFlow model, the CLM model has a simplified representation of describing the soil moisture variability while neglecting lateral water flow. The CLM model tends to overestimate the soil moisture and provided similar soil moisture estimation in grid cells that have the same soil type and plant functional type. The lateral subsurface process in CLM-ParFlow lead to soil water redistribution and improvements in prediction. The coupled model can describe the spatial variability of soil moisture. It is worth to consider lateral subsurface flow in LSMs to have more accurate soil moisture simulation.
However, some limitations should be noted. First, lateral subsurface flow takes mainly place in the saturated subsurface. However, there is also evidence for lateral flow in the top layer of the unsaturated zone for sloping soil due to rainfall dynamics. Nevertheless, the CRNS measurements provide in the first place only information on soil moisture, and are less suited to evaluate how well the influence of groundwater is represented by models. Second, the CRNS measurements might be slightly biased, caused by the limited number of observation sites, scale mismatches, and imperfect calibration. The bias is simply set-aside when using the same statistical evaluation methods in order to compare these results with other remote sensing and land surface modeling studies. Finally, this study only covers three seasons for the years of 2017 (quite average conditions) and 2018 (very dry). A longer time may be desirable to better evaluate the relative performance of the model, including different weather conditions. Also, a finer soil map resolution and larger study domain would be desirable in future studies.

Conflicts of Interest:
The authors declare no conflict of interest.