Multi-scale analysis of bias correction of soil moisture

Remote sensing, in situ networks and models are now providing unprecedented information for environmental monitoring. To conjunctively use multi-source data nom inally representing an identical variable, one must resolve biases existing between these disparate sources, and the char acteristics of the biases can be non-trivial due to spatiote mporal variability of the target variable, inter-sensor diffe rences with variable measurement supports. One such example is of soil moisture (SM) monitoring. Triple collocation (TC) bas ed bias correction is a powerful statistical method that incre asingly being used to address this issue but is only applicable to the linear regime, whereas nonlinear method of statistical moment matching is susceptible to unintended bias es originating from measurement error. Since different physi cal processes that influence SM dynamics may be distinguishable by their characteristic spatiotemporal scales, we pro pose a multi-time-scale linear bias model in the framework of a wavelet-based multi-resolution analysis (MRA). The joint MRA-TC analysis was applied to demonstrate scaledependent biases between in situ, remotely-sensed and modelled SM, the influence of various prospective bias correction schemes on these biases, and lastly to enable multi-sca le bias correction and data adaptive, nonlinear de-noising vi a wavelet thresholding.


Introduction
Global environmental monitoring requires geophysical measurements from a variety of sources and sensors to close the information gap. However, different direct and remote sensing, and model simulation can yield different estimates due to different measurement supports and errors. Soil moisture (SM) is one such variable that has garnered increasing interest due to its influences on atmospheric, hydrologic, geo-morphic and ecological processes (Rodriguez-Iturbe, 2000;GLACE Team et al., 2004;Legates et al., 2011). It also represents an archetype of the aforementioned problem, where in situ networks, remote sensing and models jointly provide extensive SM information.
In situ networks usually provide point-scale measurements; satellite retrieval of shallow SM at a mesoscale footprint of 10-50 km must resort to a homogeneity or dominantfeature assumption, whereas modelled SM depends on the simplified model parameterization, and the quality, resolution and availability of forcing data. Subsequently, the spatial (lateral and vertical) variability of SM can lead to systematically different measurements regarded as biases. Descriptive or predictive spatial SM statistics can be used to relate point-scale to mesoscale estimates (Western et al., 2002), but in situ data are often limited in describing the spatial heterogeneity of SM. However, without bias correction, it is not possible to conduct meaningful comparisons between in situ, satellite-retrieved and modelled SM for validation  and optimal data assimilation (Yilmaz and Crow, 2013). Standard bias correction methods are now increasingly being applied to SM assimilation in land models (Reichle et al., 2007;Kumar et al., 2012;Draper et al., 2012), numerical weather prediction (Drusch et al., 2005;Scipal et al., 2008a) and hydrologic models (Brocca et al., 2012).  proposed matching statistical moments of the data, while linear methods based on simple regression and matching dynamic ranges have also been considered (e.g. Su et al., 2013a). But these methods can induce artificial biases in the signal component of the corrected data as the error statistics were ignored; this also suggests a connection that the issue of bias correction is inseparable from that of error characterisation (Su et al., 2014a).
Triple collocation (TC) (Stoffelen, 1998), which is a form of instrument-variable regression (Wright, 1928; C.-H. Su and D. Ryu: Multi-scale analysis of bias correction of soil moisture 2014a), is increasingly being used to address these issues in oceanography (Caires and Sterl, 2003;Janssen et al., 2007) and hydrometeorology (Scipal et al., 2008b;Roebeling et al., 2013). In particular, it was used to estimate spatial point-tofootprint sampling errors (Miralles et al., 2011;Gruber et al., 2013), and correct biases in SM (Yilmaz and Crow, 2013). Based on an affine signal model and additive orthogonal error model, it assumes that representativity differences are manifested as additive and multiplicative biases. But these assumptions may have limited validity, as the temporal behaviour of SM may vary across different spatial scales, driven by a continuum of localised and mesoscale influences (e.g. Entin et al., 2000;Mittelbach and Seneviratne, 2012). Specifically, the coupling of SM with precipitation and evaporative losses (controlled by temperature, humidity, wind speed) varies across spatial scales. This can be more pronounced at places where surface hydrological features (e.g. topography, infiltration rate and storage capacity) are highly heterogeneous. Thus, the biases are likely to be non-systematic across short and long timescales on different spatial scales and errors are non-white, undermining the utility of the affine model. One possible remedy is to apply bias correction, either TC or statistical-moment matching, only to anomaly time series (Miralles et al., 2011;Liu et al., 2012;Su et al., 2014a), but it remains unclear how these methods affect the signal and noise components in the corrected data. Alternatively a moving time window can be used to examine the time-varying statistics of time series (Loew and Schlenz , 2011;Zwieback et al., 2013;Su et al., 2014a).
Given the possible (time)scale dependency in biases and errors, we propose an extension to TC analyses to include wavelet-based multi-resolution analysis (MRA) (Mallat, 1989) as a framework to (1) provide a fuller description of the temporal scale-by-scale relationships between coincident data sets; (2) study the influence of various prospective bias correction schemes; and (3) achieve multi-scale bias correction. To avoid excessive changes in the noise characteristics upon correction, TC can be further combined with the wavelet thresholding (Donoho and Johnstone, 1994) to (4) achieve non-linear, data-adaptive de-noising, with contrast to existing linear schemes (Su et al., 2013b). The techniques were applied to SM data from an in situ probe, satellite radiometry and land-surface model, but the proposed methods are generally enough to be applied to other geophysical variables.
The paper is organised as follows. Section 2 presents the study area over Australia and the SM data sets used in our pilot studies. Section 3 explains the theoretics behind MRA and applies it to SM, following by examination of scaleby-scale statistics in Sect. 4. Section 5 presents a new joint MRA-TC analysis framework, which is then applied to examine the influence of different bias correction schemes in Sect. 6. Importantly, both Sects. 4 and 6, using wavelet correlation, wavelet variance and scale-level TC analyses, provide evidence to support the need to extend traditional bulk and anomaly based analyses. Section 7 demonstrates the use of wavelet thresholding to de-noise satellite SM. Section 8 offers our concluding remarks.

Study areas and data sets
We consider in situ, satellite-retrieved and modelled SM over Australia. For an in-depth study, we consider pointscale and pixel-scale SM estimates at K1 monitoring site (147.56 • longitude, −35.49 • latitude) situated at Kyeamba Creek catchment, southeastern Australia (Smith et al., 2012;Su et al., 2013a). The in situ SM (INS as shorthand) was sampled at 30 min intervals, 0-8 cm depth using a timedomain interferometer-based Campbell Scientific 615 probe during November 2001-April 2011. The region experiences a temperate (Cfb) climate characterised by seasonally uniform rainfall but variable evapotranspiration forcing, so that SM varies between dry in summer (December-February) to wet in winter (June-August). The creek is located on gentle slopes with rain-fed cropping and pasture, and the soil varies from sandy to loam. Figure 1 illustrates the land cover, elevation, monthly rainfall accumulation (from 2002 to 2011), and clay content over the region.
The satellite SM was retrieved by AMSR-E (Advanced Microwave Scanning Radiometer for Earth Observing System; AMS) of the AQUA satellite. The retrieval is based on an inversion of the forward radiative transfer model of a vegetation-masked soil surface, relating observed brightness temperature to soil dielectric constant estimates. A dielectric mixing model is then used to related the dielectric constant to volumetric SM. The combined C/X-band 1/4 • × 1/4 • gridded, half-daily (∼ 1.30 a.m./p.m. LT -local time) version 5 product (July 2002-October 2011) is based on the Land Parameter Retrieval Model (Owe et al., 2008). C-band (X-band) has a shallow sampling depth of ∼ 1-2 cm (∼ 5 mm), although it is mostly C-band data over Australia due to relatively small radio frequency interference. Given the 1-2 day revisit times of the satellite, there is a significant number of missing values in the AMS data. However, we found that (not shown) over 99 % (95 %) of the gaps over Australia are ≤ 1.5 day (≤ 1 day) long. For use in wavelet analysis (Sect. 4), a one-dimensional (1-D in time) interpolation algorithm (Garcia, 2010) based on discrete cosine transform (Wang et al., 2012) was applied to infill gaps of lengths ≤ 5 days in AMSR-E. Other interpolation methods were trialled; e.g. linear interpolated AMSR-E shows great similarities to the DCT interpolated data, while cubic spline interpolation leads to spurious peaks.
The modelled SM is taken from MERRA (Modern Era Retrospective-analysis for Research and Applications) -Land produced by the Catchment land surface model GEOS version 5.7.2. The MERRA atmospheric re-analysis is driven by a vast collection of in situ observations of atmospheric and surface winds, temperature, and humidity, and remote sens- ing of precipitation and radiation (Rienecker et al., 2011). The MERRA land-only fields were post-processed by reintegrating a revised Catchment model with more realistic precipitation forcing to produce the MERRA-Land (MER as shorthand) data set . The resultant SM field corresponds to the hourly averages of the uppermost layer (0-2 cm) and is gridded on a 2/3 • × 1/2 • grid. The three data are co-located spatially via nearest neighbour and temporally at around the satellite overpass times of 1.30 a.m./p.m. LT. Their time series are plotted in blue in the first panels of Fig. 2. While co-located, the three methods observed SM dynamics over different locations and areas of the catchment (Fig. 1), due to differences in their pixel resolutions and alignments.
Continental-scale AMS and MER data over Australia are also considered. The continent has great variability in climatic and land surface characteristics. Most of the northern regions experience a tropical savannah (Aw) Köppen-Geiger climate as classified by Peel et al. (2007), central Australia is largely arid desert (BWh), and eastern mountainous areas have a temperate climate with no dry seasons (Cf). The southwestern regions similarly have a temperate climate, but with dry summers (Cs). These temperate regions have higher vegetation compared to the tropical north with moderate vegetation cover.

Multi-scale decomposition of soil moisture
The observed Kyeamba SM (denoted by blue curves p in Fig. 2) exhibits a long-term cycle of wet and dry years due to the El Niño-Southern Oscillation and seasonal and diurnal cycles originating from the fluctuations in vegetation and solar radiation, and experiences transient decay from various loss mechanisms, and abrupt increase from individual rainfall events. Their influences on observed SM can vary with the measurement methods. To unravel these differences, we turn to wavelets as the analysing kernels to study variability on individual broad-to-fine timescales. The scale under investigation is temporal for the rest of the paper, unless stated otherwise.
The 1-D orthogonal discrete wavelet transform (DWT) enables MRA of a time series p(t) of dyadic length N = 2 J and a regular sampling interval t by providing the mechanism to go from one resolution to another via a recursive function with an expectation value E(p under j 0 levels of decomposition. Loosely speaking, for a half-daily time series, the detail time series p j for j = 1, 2, 3, . . . corresponds to (fine-scale) dynamics observed on 1 day (1 d), 2 d, 4 d, etc., timescales, while the approximation time series p (a) j for j = 1, 2, 3, . . . contains (broad-scale) dynamics on scales longer than 1 d, 2 d, 4 d, etc.
In Eq.
(3), each of these components is further decomposed into a linear summation of n j = N/2 j number of basis functions φ j k and ψ j k with scale of variability 2 j t and temporal location k 2 j t. The weighting or wavelet coefficients, determined via DWT of p, measure the similarity between p and the bases via the inner products p (a) j k ≡ p, φ j k and p j k ≡ p, ψ j k . Hence the coefficients indicate changes on a particular scale and location, and enable the above scale-by-scale decomposition. Note that the bases are defined in L 2 (R) space and satisfy orthonormality conditions prescribed by φ j k , φ j k = δ j j δ kk , ψ j k , ψ j k = δ j j δ kk , φ j k , ψ j k = 0, where δ is the Kronecker delta function. For detailed expositions of the mathematical theory of wavelets and MRA, consult Daubechies (1992) and Mallat (1989).
The detail and approximated time series of Kyeamba's SM are illustrated in subsequent panels of Fig. 2, analysed using the Daubechies D(4) wavelet for j 0 = 8. On the finest scales j = 1-2 (1-2 d), the details show variability due to rainfall wetting, and over the next set of scales j = 2-5 (2-16 d) they describe transient moisture loss. The p (a) 6 (≥ 32 d) component accounts for several scales of fluctuations over seasonal, inter-annual, and long-term timescales. For comparison, the standard monthly average analyses of the original time series p are superimposed on p (a) 6 (red dots). The differences between the details of the three SM are apparent on the finest scales, with AMS and MER showing greater variability and amplitude compared to INS. However, the similarity of their temporal patterns, in both details and approximations, grows with increasing scales j > 3 (see also

Multi-scale statistics
MRA enables direct comparisons between any two representations p = {X, Y } of a given variable f (e.g. SM) on various temporal scales independently, owing to the orthonormal properties of wavelet bases. It also offers an additional degree of freedom in temporal positions (using the index k) to allow better representation of local variability. By subsetting the wavelet coefficients over certain range of k values, nonstationary statistics can also be examined. However, in this work, we consider only variability across j and assume stationarity on each scale. Pearson's linear correlation R and variance analyses (see Appendix A) are performed on the Kyeamba's INS, AMS and MER SM (as p in Eq. 2) detail (p j ) and approximation (p (a) j ) time series in Fig. 3. The strength of MRA is that since the detail time series p j on a given scale j does not contain variations on timescales greater than j , the weak-sense stationarity conditions can be better met.
Before proceeding, we recall that weak R indicate the presence of noise and/or the presence of non-linear correlation between any pairs of the data, while differences in standard deviation can also indicate the presence of noise, but also an extraneous signal and/or multiplicative bias. Typically one invokes a linearity assumption and assumes an affine relation between the signal components of the different data and an additive noise model (more later in Sect. 5), so that the differences between the data are attributed only to an overall additive bias E(X) − E(Y ), multiplicative biases, and noise. While we adopt this simplistic viewpoint here, its limitations to properly account for variable lateral and vertical measurement supports should be noted. For instance, short-timescale SM dynamics show increasing attenuation in amplitude, but are also delayed in time in deeper soil columns (e.g. Steelman et al., 2012). Additionally, SM is physically bounded between field capacity and residual content and these thresholds can vary with soil texture, location and depths. These effects can give rise to temporal autocorrelation in errors and undermine the linearity assumption be-tween coincident measures. Finally, the non-stationary characteristic of noise in satellite SM (Loew and Schlenz , 2011;Zwieback et al., 2013;Su et al., 2014a) due to e.g. dynamical land surface characteristics such as soil moisture (Su et al., 2014b), is not treated here.
With these considerations, we first examine the correlations between the three data. For the detail time series ( Fig. 3a), their correlations are lowest on the finest scales (R < 0.2) but generally improves with scale (R > 0.5), as noted previously. There is however no data pair that shows consistently higher R than other pairs: is highest on other scales. Comparing their approximation time series (Fig. 3b), R between AMS and MER are higher than the other two pairs, ranging from (j = 2) 0.8 to 0.92 (j = 8), largely due to the strong correlation between their respective p 8 and p (a) 8 . In other words: on one hand, AMS and MER both show skill in representing some aspects of the in situ SM temporal variability; on the other hand, stronger AMS-MER correlations on the coarsest (temporal) scales and their mesoscale spatial resolutions would indicate lesser representativeness of in situ measurement on these spatio-temporal scales.
Furthermore, we observe that R(p . Aside from including more noise in the approximation time series, adding components with different multiplicative biases (more later in Sect. 6) can also diminish the correlations. The scale dependence of multiplicative biases and added noise can contribute to the contrasting results of applying TC to raw versus anomaly SM time series in Draper et al. (2013). In particular, given the presence of noise in p j for j ≥ 7, error analysis of the anomaly SM (i.e. in p j for j ≤ 6) will under-estimate the total error in the raw data p.
Next, Fig. 3c plots their wavelet spectra that decompose total variance var(p) into individual scales var(p j ) ≡ std(p j ) 2 . The three data show clear differences in their standard deviation (SD) profile, both in the fine and coarse scales. As already noted, both noise and/or multiplicative biases are possible contributing factors such that noise can inflate the variance, while biases can cause suppression or inflation. Following the visual inspection of Fig. 2 and the noted weak correlations R(INS j , AMS j ) and R(INS j , MER j ) at small j , it can be argued that there is significant noise in AMS (for j = 1-3) and MER (j = 1). This in turn leads to their larger SD cf. INS. On coarser scales where R values are significantly higher, the differences in SD may be attributed more to multiplicative biases. For instance for their p 8 and p (a) 8 components, AMS and MER shows larger SD and thus positively biased relative to INS. Figure 4 extends the variance and correlation analyses between AMS and MER to the Australian continent using their coincident data from the period July 2002-October 2011. The spatial maps of SD differences ( SD) and correlations show significant variability in the statistics with timescales and spatial locations. On the finest scale j = 1, the similarity between the difference map (Fig. 4a) and the TC-derived error map of AMSR-E (see Fig. 6a in Su et al., 2014a) in terms of spatial variability and the low AMS-MER correlations (Fig. 4f) support our observation that the detail time series AMS 1 is noise dominated. Weak negative correlation between AMS 1 and MER 1 can also be observed over arid regions. By contrast, owing to the strong correlation R ∼ 0.6-0.9 ( Fig. 4g and h) on the coarse scales, the causes of SD ( Fig. 4c and d) are related to biases. In particular, at j > 8, the SD map in Fig. 4d also suggests a possible association between biases and climatology or land cover characteristics, with negative biases dominating northern tropical (Aw) and semi-arid (BS) regions, and positive biases in temperate, vegetated regions (Cs and Cf) over southeastern and southwest-ern Australia. The visual comparisons between scale-level SD and bulk SD enable stratification of the continent to central arid regions of higher noise identified in j = 1 and 2 and temperate (tropical) regions, with a positive (negative) bias seen on coarser scales.

Joint MRA-TC analysis
In order to quantify observed differences between the data, we propose a scale-dependent linear model: a multi-scale (MS) model that distinguishes the signal components of the two data X and Y via an overall additive bias and a set of positive scaling coefficients α p,j , α p , and assumes an additive and zero-mean independent but non-white noise model p (t). Focusing on the zero-mean signal and noise components, the "structural relationship" model reads where the signal and noise components have been decomposed into their multi-resolution forms. The standard assumptions of orthogonal and mutually uncorrelated errors are used, so that the covariance cov(f j , e p,j ) = 0, cov(f , e p ) = 0, cov(e p,j , e q,j ) = 0, cov(e p,j , e q ) = 0 and cov(e p , e q ) = 0 for p = q, p, q,∈ {X, Y }. The differences in the values of the scaling coefficients between data, i.e. α X,j = α Y,j , signify multiplicative biases on individual scales. To see this, we express their mean-squared deviation MSD ≡ E[(Y − X) 2 ] in terms of variables in Eqs. (4) and (5) to arrive at The first term is the additive bias, and the summation consists of scale-specific multiplicative biases proportional to (α X,j − α Y,j ) 2 and noise contributions from each datum. The interpretation of the discrepancies between X and Y can vary depending on the time period of the data and the analysis, and the adopted signal/noise model. By using the entire 9 year record of INS, AMS and MER data in MRA, the MS model does not observe a time-varying additive bias (e.g. from using the moving-window approach of Su et al., 2014a) or autocorrelated errors (from using the lagged covariance in Zwieback et al., 2013). Rather, MRA and the MS model enable a description of the systematic differences based wholly in terms of multiplicative biases at individual timescales, and the random differences in terms of additive noise. Specifically, this contrasts with the short time-window approach (e.g. ≤ 32 d), where multiplicative biases existing on coarse scales (p (a) 6 ) will manifest as both time-varying additive and multiplicative biases. Importantly, the model allows for different scaling coefficients between scales, i.e. α p,j = α p,j for j = j , as a form of non-linearity with f . The equality α p,j = α p = α p is therefore a special case of (bulk) linearity. As our focus of the above model is the multiplicative biases and noise, for convenience of notation, we remove the mean of the X and Y prior to MRA and bias correction. Furthermore, without the loss of generality, we choose X as the reference henceforth and let α X,j , α X = 1.
By using a third independently derived representation (Z) of f , TC enables estimation of the required scaling coefficients and noise std( p,j ) (Appendix B). As we will see later, these estimates are needed for bias correction and de-noising. Within the operating assumptions of TC, TC estimates are unbiased and consistent; that is, the estimatedα Y,j = α Y,j as the asymptotic limit. However, TC's superiority is dependent on the availability of a strong instrument and a large sample for statistical analyses (Zwieback et al., 2012;Su et al., 2014a). Standard linear estimators, namely ordinary leastsquare (OLS) regression and variance matching (VAR), can be considered as substitutes, although they are biased estimators of α when X and Y are both noisy (Yilmaz and Crow, 2013;Su et al., 2014a), e.g. OLS yieldsα Y,j < α Y,j . In summary, we propose that combining these estimators with MRA via the MS linear model enables investigation into the distribution of the multiplicative biases and additive noise over j , and their response to various bias correction schemes.

Multi-scale analysis of bias correction
Consider now the bias correction of Y to produce a corrected datum Y * that "matches" X. Different interpretations of a "match" and assumptions about signal and noise statistics lead to different bias correction schemes. To describe matching, there are different choices of optimality criterion. The first is based on matching the statistics of the signalonly component of Y * to that of X. This approach requires consistent estimation of slope parameters α's and the resultant statistics of X and Y * may differ due to different noise statistics. The second is based on the matching of the statistical moments between Y * and X (e.g. VAR matching), although the statistics of their constitutive signal components may differ for the same reason. The third is based on the minimum-variance principle of minimizing the least-square difference between Y * and X (i.e. the OLS estimation), but as already noted, the estimator becomes inconsistent when there are measurement errors in X and Y . Following our theoretical model in Sect. 5, we define our optimality criterion based on the first criterion of matching the first two moments of the signal components in X and Y so that Y * is suitable for bias-free data assimilation. In particular, Yilmaz and Crow (2013) have shown that residual multiplicative biases due to a sub-optimal bias correction scheme will cause filter innovations to contain residual signal and sub-optimal filter performance. Thus, within the paradigm of the MS model, our goal of bias correction is to minimise the difference |α Y * ,j − 1| for α X,j = 1, so that the multiplicative bias terms in Eq. (6) are eliminated.
-Bulk linear rescaling assumes bulk linearity between X and Y so that the correction equation is whereα Y is given by TC for our objective. When the bulk linearity is satisfied, this approach ensures that the statistical properties (SD and higher moments) of the signal components in X and Y * are identical. Linear rescaling usingα Y values estimated by OLS and VAR matching have previously been considered by e.g. Su et al. (2013a); but due to error-in-variable biases, they can induce artificial biases in the signal component of Y * even if the bulk linearity condition is valid.
-Bulk cumulative distribution function (CDF) matching assumes non-linearity between X and Y and transforms Y * so that  cdf Y * = cdf(X), where cdf(•) computes the CDF. This ensures that the mean, SD, and higher statistical moments of X and Y * are identical, but the statistical properties of their signal and noise components that make up X and Y * are 24 C.-H. Su and D. Ryu: Multi-scale analysis of bias correction of soil moisture not necessarily identical. In particular, when the relative signal and noise statistics in the two data are different, CDF matching leads to artificial biases between the signal components in X and Y * . As with VAR matching of first two moments, the CDF counterpart is expected to contain extraneous contribution of the noise variances in the mapping of the second moment, as well as at higher moments (Su et al., 2014a). The issue can be exacerbated by variable signal and noise statistics on different scales.
-Anomaly/seasonal (A/S) linear rescaling allows biases between X and Y to be different on two scales of variation. In practice, the useful information content in observations is primarily based on their representation of anomalies, where observations are assumed in a particular land surface model's unique climatology (Koster et al., 2009). The correction is therefore limited to the anomalies, although other components (e.g. seasonal fluctuation and long-term trend) may be preserved to validate model prediction. Here the linear correction using TC estimator is applied to match the characteristics of each component -anomaly (i = A) and seasonal (S) -separately, so that the corrected Y has the form In one approach, p S is computed using moving window averaging of multiyear data within a window size of 31 days centered on a given day of year (Miralles et al., 2011;Su et al., 2014a), so that inter-annual cycles and long-term trends are retained in p A . In an alternative approach (Albergel et al., 2012), a sliding 31 day window is used such that p j for half-daily time series. In this work, the former, more conventional approach was taken.
-A/S CDF matching applies CDF matching to anomaly and seasonal components separately as per Eq. (9) but with cdf(Y * i ) = cdf(X i ). The application of CDF matching to the anomaly component of soil moisture data was considered by Liu et al. (2012).
In relation to Eq. (6), this approach obviously eliminates that the multiplicative terms in the summation. The bulk and A/S linear correction schemes can be considered as special cases of MS rescaling where information from multiple scales are aggregated and corrected jointly. Other aggregations of the information from different subsets of scales are also possible, but they will similarly be conceived based on one's understanding or assumptions of the underlying specific processes driving SM dynamics. Investigations into suitable aggregations are beyond the scope of this work, hence we implemented the most elaborate decomposition. If joint linearity exists between two or more scales, their α Y,j values will be similarly valued for use in Eq. (10).
For illustrations, we correct the biases in AMS and MER SM with respect to INS SM at Kyeamba using the above five schemes. Using the above notations, AMS and MER are treated as Y , the corrected AMS * and MER * as Y * , and INS as X. MRA-TC was applied to observe their consequences in Fig. 5. In the upper panel, estimatedα Y,j andα Y * ,j values provide diagnostics for detecting the presence of multiplicative biases before and after application of the correction schemes. The lower panel plots the SD of Y j and Y * j and their associated noises Y,j and Y * ,j . The values of the scaling coefficients α Y,j (before correction) and α Y * ,j (after), and the noise std( Y,j ) and std( Y * ,j ) were estimated using TC. But where TC estimates could not be retrieved (for j = 1-2) due to negative correlation amongst the data triplet (e.g. resulting from significant noise and weak instrument), OLS-derived (under) estimates serve as a guide for the above diagnostic purposes. Similarly, the total SD is a guide for noise SD in these cases. Figure 5a shows the MRA of the biases and noise in the pre-corrected data Y . There is considerable variability inα Y,j across the scales, ranging from 0.5 to 1.8 for AMS, and from 0.5 to 1.4 for MER. In particular, theirα Y andα Y,8 deviate significantly from 1, and are responsible for the larger SD (cf. INS) observed in Fig. 3c. Biases also exist on almost all other scales of AMS and MER. In the lower panel, the values of std( Y,j ) relative to std(Y j ) indicate the dominance of noise in the small scales j = 1-3. This explains the low R values between AMS (and MER) and INS in Fig. 3a. Furthermore, the signal-to-noise ratios are variable with scales and data sets, highlighting the importance of using a correction scheme that takes the signal-vs.-noise statistics into considerations. The TC-based scheme is limited to the linear case, and the CDF scheme ignores such a variability.
The MRA of the corrected data Y * are shown in Fig. 5b-f. In addition we assess the level of agreement between corrected AMS * and INS time series in Table 1 using their rootmean-squared deviation (RMSD) and correlation R. The time series plots are shown in Fig. 6 to support interpretations. These additional results focus on the AMS-INS pair that best illustrates the influence of noise in AMS.
The results of bulk, A/S and MS linear rescaling can be readily interpreted. For bulk (Fig. 5b) and A/S linear (Fig. 5d) rescaling, the values ofα Y andα Y i used for their implementation (Eqs. 7 and 9) are listed in the figure. As these values are greater than unity for both AMS and MER, this leads to the Hydrol. Earth Syst. Sci., 19, 17-31  suppression of the associated signal, as well as noise, components: std(Y * j ) < std(Y j ), and std( Y * ,j ) < std( Y,j ). For AMS, the bulk linear scheme corrects the coarse-scale bias in Y (a) 8 component and rescales the noise variance, reducing RMSD from 0.09 to 0.06 m 3 m −3 . However, the fine-scale biases in Y * j are still present, and increased on some scales, e.g. at j = 4, 7 for AMS * . Additionally for A/S linear rescaling, R(AMS * ,INS) value does not change significantly and the noise are still clearly visible in Fig. 6b and d. By construction, the MS rescaling uses the estimatedα Y,j values from Fig. 5a to correct bias on all the scales. Fig. 5f shows the analysis of MS-corrected Y * . The equivalencê α Y * ,j = 1 indicates that the multiplicative biases are eliminated at j > 2. At j = 1-2, as the scaling coefficients cannot be estimated by TC, CDF matching was applied to these scales such that the biases are still present on these scales. Amid the reduction of biases, we also observed noise amplification (i.e. std( Y * ,j ) >std( Y,j )) in AMS * at j = 3, 7 and in MER * at j = 3-7, because of rescaling with less-than-unitŷ α Y,j values in Eq. (10). Indeed it is evident from Eq. (6) that it is possible to increase the noise variance and MSE when reducing the bias component of the MSE. This in turn leads to larger disagreement between INS and AMS in terms of RMSD and R, and the increased amplitudes of the noise observed in AMS in Fig. 6f.
The bulk and A/S CDF methods produced very similar results with each other, and also with their linear counterparts. There is signal and noise suppression, but the scale-level biases are retained. The signal components of Y * are negatively biased at j = 3-7 and positively biased at j = 8. The CDFcorrected AMS * shows slightly better RMSD and R with INS, owing to the reduced noise variance and a reduced bias at AMS (a) *

.
In summary, the MRA of the bulk and A/S schemes highlights the deficiency of using a correction scheme that does not take into account the scale variability of bias and the differences in noise statistics between the two data. The im- provements in RMSD and correlation between the corrected Y * and the reference X are somewhat superficial, masking the fact that the bias correction is limited to the coarsest scales. On the other hand, the A/S-based and MS methods can modify the original noise profiles in the data across the scales, by amplifying (or suppressing) noise in individual components (either Y j , Y S , or Y A ) with less-than (greaterthan) unity pre-correction α. This may be considered undesirable for an objective to produce more physically representative data with a simple error structure on the whole. Therefore, arguably, none of these methods is entirely satisfactory, in manners of not removing the multiplicative biases completely and/or changing error characteristics. From this viewpoint, the task of bias correction is seen as inseparable from that of noise reduction when considering MS (or A/S) bias correction, unless certain components in MRA were explicitly ignored.

Combining bias correction with wavelet de-noising
The last example presents an impetus to consider noise removal prior to bias correction and produce a simpler error structure in the bias-corrected data Y * . Critically, TC provides noise and signal estimates that can be used for denoising through thresholding of wavelet coefficients p j k . The basic rationale for wavelet thresholding (WT) is that insignificant detail coefficients are likely due to noise, while significant ones are related to the signal component. Thus, a coefficient is eliminated if its magnitude is less than a given threshold λ p ; otherwise, it is modified according to a transformation function (p j k ) to remove the influence of the noise (Donoho and Johnstone, 1994). One commonly used transformation is soft thresholding (Donoho, 1995), where the coefficients are modified according to λ p p j k = sign p j k max |p j k − λ p |, 0 .
Such de-noising filters have near-optimal properties in the minmax sense. We follow the BayesShrink rule of Chang et al. (2000) to define a set of scale-dependent threshold values using where the variances are provided by TC (Appendix B). This choice of threshold is near optimal under the assumption that the signal is generalised Gaussian distributed and the noise is Gaussian. When the threshold value for j = 1-2 could not be estimated using TC, CDF matching was applied. While TC is an ideal error estimator, alternative estimators for the threshold values are also available to make the de-noising a stand-alone process (Donoho and Johnstone, 1994;Donoho, 1995). After WT, the de-noised time series is constructed via inverse DWT of the modified coefficients, and can be subsequently corrected for biases. Combining with the MS bias correction scheme, biased-corrected, de-noised data are generated via The prescription, which is essentially a two-stage operation, was applied to AMS for comparisons with the previous results. The first stage of de-noising leads to smoothing of the time series, improved R with INS by 0.05, and reduced RMSD by 0.02 m 3 m −3 . The actual SM variability has become more apparent in Fig. 6g. Over-smoothing can occur due to our inability to properly distinguish signal from noise in AMS 1 and AMS 2 where the signal-to-noise ratio is very low. However, without the second stage of bias correction, the dynamic ranges of de-noised AMS and INS are visibly different, such that the improvement in RMSD with INS is limited. Combining WT and MS leads to improvement in both metrics of RMSD = 0.048 m 3 m −3 and R = 0.711, with Fig. 6h confirming that the reduced noise was not amplified by the MS rescaling.

Conclusions
This work combines MRA and TC in a new analysis framework with increased capacity to provide a more comprehensive view of the inter-data relations on short and long timescales. TC (or CDF) rescaling can be exploited on individual scales to reduce scale-specific multiplicative biases, and provide "prior" knowledge of noise for calibrating a WTbased de-noising filter. As a demonstration of principle, these methods are applied to SM data from in situ and satellite sensors and a land surface model. Using MRA, we found that the three data exhibit significantly different wavelet spectra and variable degrees of agreement on different timescales. On fine scales, the contribution of noise is most prominent, undermining the correlation between the data sets. By contrast, the biases are most apparent on coarse scales. Furthermore, these biases are non-systematic across timescales in the study region and across spatial locations over Australia, and the signal-to-noise ratios vary with scales and between the various data, pointing to the need to use correction schemes that are capable of handling such complexities.
These observations raised concerns about the possible inadequate treatment of SM data in the linear regime, even with anomaly/seasonal decomposition. Scale-by-scale linear rescaling based on a MRA-TC analysis framework offers a more comprehensive treatment of different biases on different scales, but error characteristics are found to be modified by variable rescaling, and can lead to undesirable noise amplification. The method of removing biases and noise on individual scales offers a remedy, although a few caveats should be noted. First, TC analysis requires a strong instrument and large sample, and in cases where these prerequisites are not met, we resort to sub-optimal estimation and rescaling methods. Second, the issue of non-stationarity in errors and scaling has not been addressed so far, and this can lead to biased estimates of the correction parameters for rescaling and de-noising. Despite this, DWT offers additional degree of freedom in translation parameter k to accommodate nonstationarity. Third, given the theoretic viewpoint presented in this work, further evaluations based on assimilation of data treated by different schemes are still warranted to assess their practical impacts. Notwithstanding these factors, MRA-TC analysis can be an important tool to allow better characterisation of the inter-sensor differences and to develop more effective strategies in harmonising a broad range of observational data records in oceanography and hydrometeorology. MRA enables the (bulk) variance var(p) of a time series p to be decomposed into wavelet variances var(p j ) on different scales j . Analogous to a Fourier spectrum, the expansion of var(p) yields a wavelet spectrum and is given by where the variance of the approximation time series p (a) j 0 can be expressed in terms of that of the detail time series p j .
Similarly, wavelet covariance cov(X j , Y j ) at a given j indicates the contribution of covariance between two time series (X, Y ) on that scale. Specifically, the wavelet covariance on scale j can be expressed as cov X j , Y j = 1 n j n j k=1 X j k Y j k , noting that there is an equivalence of computing (co)variance in the wavelet and time domains. To exclude the boundary influence of a finite-length time series and missing values in the time series, an estimator of the wavelet covariance can be constructed by excluding the coefficients affected by the boundaries and gaps, followed by renormalisation. In the paper, we find it more intuitive to report the wavelet correlation, namely R X j , Y j = cov X j , Y j var X j var Y j . (A4)