Pilot-point based multi-objective calibration in a surface–subsurface distributed hydrological model

Abstract

In this study, a fully-coupled surface–subsurface, distributed, physics-based hydrological model was calibrated using the pilot-point method. A minimum variance field rule was included in the objective function to regularize the extensive calibration exercise that included 74 parameters (72 associated with pilot points and two spatially-invariant channel parameters). Because the overland and vadose zone systems are not in permanent hydrological connection, the information contained in the observation points may not be accessible by the pilot points at all times, rendering them insensitive to the observations and hindering the calibration process. An analysis of the spatial and temporal variability of parameter sensitivities was done to explore how the information contained in local observations spreads from the observation points to the pilot points, where parameter values are identified. The results show that the channel flow time series is valuable to identify the parameters at all pilot-point locations, indicating that the information in channel flow propagates to the entire basin. However, information in soil moisture measurements is of local extent and thus only valuable to identify the parameters at locations close to the observation point.

Editor D. Koutsoyiannis; Associate editor I. Nalbantis

Citation Maneta, M.P. and Wallender, W.W., 2013. Pilot-point based multi-objective calibration in a surface–subsurface distributed hydrological model. Hydrological Sciences Journal, 58 (2), 390–407.     

An advanced process for evaluating a linear dielectric constant–bulk electrical conductivity model using a capacitance sensor in field conditions

Abstract

The Hilhorst model was used to convert bulk electrical conductivity (σ_b) to pore water electrical conductivity (σ_p) under laboratory conditions by using the linear relationship between the soil dielectric constant (ε_b) and σ_b. In the present study, applying the linear relationship ε_b–σ_b to data obtained from field capacitance sensors resulted in strong positive autocorrelations between the residuals of that regression. We were able to derive an accurate offset of the relationship ε_b–σ_b and to estimate the evolution of σ_p over time by including a stochastic component to the linear model, rearranging it to a time-varying dynamic linear model (DLM), and using Kalman filtering and smoothing. The offset proved to vary for each depth in the same soil profile. A reason for this might be the changes in soil temperature along the soil profile.

Editor D. Koutsoyiannis; Associate editor M.D. Fidelibus     

In situ lateral-loading test and micro–macro-scale simulation of an existing RC frame

An in situ pseudo-dynamic lateral loading test was conducted on a real three-storey, one-by-one bay RC frame structure to study its seismic performance. The test damage process shows that it is a column-sway collapse mechanism with the damage concentrated especially at the ends of the bottom columns. In accordance with the comparison of the macro-scale and micro-scale finite element models, an element-coupling model (ECM) is proposed to acquire the micro-scale damage details with suitable calculation cost. In the ECM, macro-scale and micro-scale elements are coupled together by the multi-point constraint equations, and a plastic-damage model with modified yield surface is used to describe the nonlinear behavior of the masonry material. The calculated results show that the ECM could correctly reflect the global behavior and the micro-scale damage details with efficient cost, and the relationship of the local damage and global performance could be obtained, which is helpful to study the damage mechanism.     

Simulating pathways of subsurface oil in the Faroe–Shetland Channel using an ocean general circulation model

Little is known about the fate of subsurface hydrocarbon plumes from deep-sea oil well blowouts and their effects on processes and communities. As deepwater drilling expands in the Faroe–Shetland Channel (FSC), oil well blowouts are a possibility, and the unusual ocean circulation of this region presents challenges to understanding possible subsurface oil pathways in the event of a spill. Here, an ocean general circulation model was used with a particle tracking algorithm to assess temporal variability of the oil-plume distribution from a deep-sea oil well blowout in the FSC. The drift of particles was first tracked for one year following release. Then, ambient model temperatures were used to simulate temperature-mediated biodegradation, truncating the trajectories of particles accordingly. Release depth of the modeled subsurface plumes affected both their direction of transport and distance travelled from their release location, and there was considerable interannual variability in transport.     

Seasonal drought predictability in Portugal using statistical–dynamical techniques

Atmospheric forecasting and predictability are important to promote adaption and mitigation measures in order to minimize drought impacts. This study estimates hybrid (statistical–dynamical) long-range forecasts of the regional drought index SPI (3-months) over homogeneous regions from mainland Portugal, based on forecasts from the UKMO operational forecasting system, with lead-times up to 6 months. ERA-Interim reanalysis data is used for the purpose of building a set of SPI predictors integrating recent past information prior to the forecast launching. Then, the advantage of combining predictors with both dynamical and statistical background in the prediction of drought conditions at different lags is evaluated. A two-step hybridization procedure is performed, in which both forecasted and observed 500 hPa geopotential height fields are subjected to a PCA in order to use forecasted PCs and persistent PCs as predictors. A second hybridization step consists on a statistical/hybrid downscaling to the regional SPI, based on regression techniques, after the pre-selection of the statistically significant predictors. The SPI forecasts and the added value of combining dynamical and statistical methods are evaluated in cross-validation mode, using the R² and binary event scores. Results are obtained for the four seasons and it was found that winter is the most predictable season, and that most of the predictive power is on the large-scale fields from past observations. The hybridization improves the downscaling based on the forecasted PCs, since they provide complementary information (though modest) beyond that of persistent PCs. These findings provide clues about the predictability of the SPI, particularly in Portugal, and may contribute to the predictability of crops yields and to some guidance on users (such as farmers) decision making process.     

Radial reactive solute transport in an aquifer–aquitard system

Radial reactive transport is investigated in an aquifer–aquitard system considering the important processes such as advection, radial and vertical dispersions for the aquifer, vertical advection and dispersion for the aquitards, and first-order biodegradation or radioactive decay. We solved the coupled governing equations of transport in the aquifer and the aquitards by honoring the continuity of concentration and mass flux across the aquifer–aquitard interfaces and recognizing the concentration variation along the aquifer thickness. This effort improved the averaged-approximation (AA) model, which dealt with radial dispersion in an aquifer–aquitard system by excluding the aquitard advection. To compare with our new solution, we expanded the AA model by including the aquitard advection. The expanded AA model considerably overestimated the mass in the upper aquitard when an upward advection existed there. The rates of mass change in the upper aquitard from the new solution and the AA model solution increased with time following sub-linear fashions. The times corresponding to the peak values of the residence time distributions for the AA model, the expanded AA model, and the new model were almost the same. The residence time distributions seemed to follow the Maxwell–Boltzmann distribution closely when plotting the time in logarithmic scale. In addition, we developed a finite-element COMSOL Multiphysics simulation of the problem, and found that the COMSOL solution agreed with the new solution well.     

On the importance of soil moisture in calibration of rainfall–runoff models: two case studies

Streamflow modelling results from the GR4H and PDM hydrological models were evaluated in two Australian sub-catchments, using (1) calibration to streamflow and (2) joint-calibration to streamflow and soil moisture. Soil moisture storage in the models was evaluated against soil moisture observations from field measurements. The PDM had the best performance in terms of both streamflow and soil moisture estimations during the calibration period, but was outperformed by GR4H during validation. It was also shown that the soil moisture estimation was improved significantly by joint-calibration for the case where streamflow and soil moisture estimations were poor. In other cases, addition of the soil moisture constraint did not degrade the results. Consequently, it is recommended that GR4H be used, in preference to the PDM, in the foothills of the Murrumbidgee catchment or other Australian catchments with semi-arid to sub-humid climate, and that soil moisture data be used in the calibration process.     

Soil–pile-bridge structure interaction in liquefying ground using shake table testing

The evaluation of seismic pile response is particularly useful for geotechnical engineers involved in the design of foundations in liquefying site. Shake table testing was performed to study the dynamic interactive behavior of soil–pile foundations in liquefying ground under different shaking frequency and amplitude. The soil profile consisted of a clayey layer over liquefiable sand over clay. The model was tested with a series of El Centro earthquake motions with peak accelerations ranging from 0.15g to 0.50g, and time step from 0.006 to 0.02 s. Representative data, including time histories of accelerations and excess pore pressure ratios that characterize the important aspects of soil–pile interaction in liquefying ground are presented. The shaking frequency has no significant effect on the magnitudes of excess pore pressure ratio, ground and pile accelerations and pile bending moments. Excess pore pressure ratio, ground acceleration and pile acceleration, and pile bending moment largely depend on the shaking amplitude.     

An in situ 14C–10Be Bayesian isochron approach for interpreting complex glacial histories

We present a Bayesian isochron approach to interpret measurements of multiple cosmogenic nuclides from glacially modified bedrock surfaces with complex exposure histories. An isochron approach explicitly incorporating glacial erosion is ideally suited for this problem; such erosion must be accounted for but has traditionally been ignored. Previous methods required treating each sample individually (to account for glacial erosion) and subsequently averaging results for the entire dataset. Geological considerations, however, suggest a more robust approach is to treat samples in the dataset here (and samples from other conceivable datasets) simultaneously. The Bayesian isochron method is applied to a previously published set of in situ ¹⁴C and ¹⁰Be measurements from a set of samples spanning the forefield of the Rhone Glacier, Switzerland. Results indicate 6.4 ± 0.5 kyr of integrated exposure and 4.7 ± 0.5 kyr of cumulative burial, similar to previous estimates, but with much smaller uncertainties. The reduced uncertainties result from fitting the exposure and burial duration to the entire dataset, while explicitly accounting for glacial erosion. The method presented here should be applicable with minor modifications in a number of geologic settings, and further demonstrates the utility of paired in situ ¹⁰Be and ¹⁴C measurements for unraveling complex exposure histories over during the Holocene and late Pleistocene.     

Simulating temporal variability in catchment response using a monthly rainfall–runoff model

Abstract

Temporal variability can result from shifts in climate, or from changes in the runoff response due to land- or water-use changes, and represents a potential source of uncertainty in calibrating hydrological models. Parameter values were determined using Monte Carlo parameter sampling methods for a monthly rainfall–runoff model (Pitman model) for different sub-periods on four catchments, with different types and degrees of temporal variability, in Australia and Africa. For some catchments, parameters were not dependent upon the sub-period used and fell within expected ranges given the relatively high degree of model equifinality. In other catchments, dependencies can be identified that are associated with signals contained within the sub-periods. While the Pitman model is relatively robust in the face of temporal variability, it is concluded that better simulations will always be obtained from calibration data that include signals representing the total variability in climate, land-use change and catchment responses.     

Time–frequency characterization of sub-divisional scale seasonal rainfall in India using the Hilbert–Huang transform

Time–frequency characterization is useful in understanding the nonlinear and non-stationary signals of the hydro-climatic time series. The traditional Fourier transform, and wavelet transform approaches have certain limitations in analyzing non-linear and non-stationary hydro-climatic series. This paper presents an effective approach based on the Hilbert–Huang transform to investigate time–frequency characteristics, and the changing patterns of sub-divisional rainfall series in India, and explored the possible association of monsoon seasonal rainfall with different global climate oscillations. The proposed approach integrates the complete ensemble empirical mode decomposition with adaptive noise algorithm and normalized Hilbert transform method for analyzing the spectral characteristics of two principal seasonal rainfall series over four meteorological subdivisions namely Assam-Meghalaya, Kerala, Orissa and Telangana subdivisions in India. The Hilbert spectral analysis revealed the dynamic nature of dominant time scales for two principal seasonal rainfall time series. From the trend analysis of instantaneous amplitudes of multiscale components called intrinsic mode functions (IMFs), it is found that both intra and inter decadal modes are responsible for the changes in seasonal rainfall series of different subdivisions and significant changes are noticed in the amplitudes of inter decadal modes of two seasonal rainfalls in the four subdivisions since 1970s. Further, the study investigated the links between monsoon rainfall with the global climate oscillations such as Quasi Bienniel Oscillation (QBO), El Nino Southern Oscillation (ENSO), Sunspot Number (SN), Atlantic Multidecadal Oscillation (AMO) etc. The study noticed that the multiscale components of rainfall series IMF1, IMF2, IMF3, IMF4 and IMF5 have similar periodic structure of QBO, ENSO, SN, tidal forcing and AMO respectively. As per the seasonal rainfall patterns is concerned, the results of the study indicated that for Assam-Meghalaya subdivision, there is a likelihood of extreme rare events at ~0.2 cycles per year, and both monsoon and pre-monsoon rainfall series have decreasing trends; for Kerala subdivision, extreme events can be expected during monsoon season with shorter periodicity (~2.5 years), and monsoon rainfall has statistically significant decreasing trend and post-monsoon rainfall has a statistically significant increasing trend; and for Orissa subdivision, there are chances of extremes rainfall events in monsoon season and a relatively stable rainfall pattern during post-monsoon period, but both monsoon and post-monsoon rainfall series showed an overall decreasing trend; for Telangana subdivision, there is a likelihood of extreme events during monsoon season with a periodicity of ~4 years, but both monsoon and post-monsoon rainfall series showed increasing trends. The results of correlation analysis of IMF components of monsoon rainfall and five climate indices indicated that the association is expressed well only for low frequency modes with similar evolution of trend components.     

Modelling of fluid–solid interactions using an adaptive mesh fluid model coupled with a combined finite–discrete element model

Fluid–structure interactions are modelled by coupling the finite element fluid/ocean model ‘Fluidity-ICOM’ with a combined finite–discrete element solid model ‘Y3D’. Because separate meshes are used for the fluids and solids, the present method is flexible in terms of discretisation schemes used for each material. Also, it can tackle multiple solids impacting on one another, without having ill-posed problems in the resolution of the fluid’s equations. Importantly, the proposed approach ensures that Newton’s third law is satisfied at the discrete level. This is done by first computing the action–reaction force on a supermesh, i.e. a function superspace of the fluid and solid meshes, and then projecting it to both meshes to use it as a source term in the fluid and solid equations. This paper demonstrates the properties of spatial conservation and accuracy of the method for a sphere immersed in a fluid, with prescribed fluid and solid velocities. While spatial conservation is shown to be independent of the mesh resolutions, accuracy requires fine resolutions in both fluid and solid meshes. It is further highlighted that unstructured meshes adapted to the solid concentration field reduce the numerical errors, in comparison with uniformly structured meshes with the same number of elements. The method is verified on flow past a falling sphere. Its potential for ocean applications is further shown through the simulation of vortex-induced vibrations of two cylinders and the flow past two flexible fibres.     

Modelling of fluid–solid interactions using an adaptive mesh fluid model coupled with a combined finite–discrete element model

Fluid–structure interactions are modelled by coupling the finite element fluid/ocean model ‘Fluidity-ICOM’ with a combined finite–discrete element solid model ‘Y3D’. Because separate meshes are used for the fluids and solids, the present method is flexible in terms of discretisation schemes used for each material. Also, it can tackle multiple solids impacting on one another, without having ill-posed problems in the resolution of the fluid’s equations. Importantly, the proposed approach ensures that Newton’s third law is satisfied at the discrete level. This is done by first computing the action–reaction force on a supermesh, i.e. a function superspace of the fluid and solid meshes, and then projecting it to both meshes to use it as a source term in the fluid and solid equations. This paper demonstrates the properties of spatial conservation and accuracy of the method for a sphere immersed in a fluid, with prescribed fluid and solid velocities. While spatial conservation is shown to be independent of the mesh resolutions, accuracy requires fine resolutions in both fluid and solid meshes. It is further highlighted that unstructured meshes adapted to the solid concentration field reduce the numerical errors, in comparison with uniformly structured meshes with the same number of elements. The method is verified on flow past a falling sphere. Its potential for ocean applications is further shown through the simulation of vortex-induced vibrations of two cylinders and the flow past two flexible fibres.

     

Assessing non-steady-state erosion processes using paired 10Be–26Al in southeastern Tibet

Quantifying erosion rates over various spatial and temporal scales across the Tibetan Plateau and its surrounding mountains is crucial to understanding the topographic evolution of the orogen. In this work, we report a new dataset of ¹⁰Be-derived basin-wide erosion rates from the main tributaries and streams draining the eastern Himalayan syntaxis. The 22 basin-wide erosion rates ranged from 78 ± 7 m Myear⁻¹ to 3,490 ± 612 m Myear⁻¹ across the study area. ²⁶Al was contemporarily measured to evaluate the impact of sediment storage and non-steady-state erosion processes in the syntaxis region. The paired study of ¹⁰Be and ²⁶Al reveals that several samples violated the steady-state erosion assumption and were compatible with the scenario of perturbation of reworked sediments or deeply sourced materials introduced by landslides. For most samples, deep-sourced materials with higher ²⁶Al/¹⁰Be ratios were no longer perturbing the ¹⁰Be signals in river sediments. It is possible that the deep-sourced materials had been wiped out of the basins before the collection of samples in this work. However, the perturbation of reworked sediments was observed over a range of basin scales, limiting the use of a single sediment sample as a representative erosion product for upstream basins. Compared with tectonically stable regions, the incorporation of reworked fluvial sediments, deeply sourced materials or sub-glacial eroded materials into sampled sediments led to the decoupling between basin-wide erosion rates and topographic or climatic indices. Caution should be taken when deriving erosion rates from rapidly eroding regions with old, deeply buried sediments such as the eastern Himalayan syntaxis, where calculated erosion rates may be highly overestimated.     

Soil moisture estimation using tomographic ground penetrating radar in a MCMC–Bayesian framework

In this study, we focus on a hydrogeological inverse problem specifically targeting monitoring soil moisture variations using tomographic ground penetrating radar (GPR) travel time data. Technical challenges exist in the inversion of GPR tomographic data for handling non-uniqueness, nonlinearity and high-dimensionality of unknowns. We have developed a new method for estimating soil moisture fields from crosshole GPR data. It uses a pilot-point method to provide a low-dimensional representation of the relative dielectric permittivity field of the soil, which is the primary object of inference: the field can be converted to soil moisture using a petrophysical model. We integrate a multi-chain Markov chain Monte Carlo (MCMC)–Bayesian inversion framework with the pilot point concept, a curved-ray GPR travel time model, and a sequential Gaussian simulation algorithm, for estimating the dielectric permittivity at pilot point locations distributed within the tomogram, as well as the corresponding geostatistical parameters (i.e., spatial correlation range). We infer the dielectric permittivity as a probability density function, thus capturing the uncertainty in the inference. The multi-chain MCMC enables addressing high-dimensional inverse problems as required in the inversion setup. The method is scalable in terms of number of chains and processors, and is useful for computationally demanding Bayesian model calibration in scientific and engineering problems. The proposed inversion approach can successfully approximate the posterior density distributions of the pilot points, and capture the true values. The computational efficiency, accuracy, and convergence behaviors of the inversion approach were also systematically evaluated, by comparing the inversion results obtained with different levels of noises in the observations, increased observational data, as well as increased number of pilot points.     

Lineaments in the Shamakhy–Gobustan and Absheron hydrocarbon containing areas using gravity data

In this study, we purposed to investigate the edge of geostructures and position of existing faults of the Shamakhy–Gobustan and Absheron hydrocarbon containing regions in Azerbaijan. For this purpose, the horizontal gradient, analytic signal, tilt angle, and hyperbolic of tilt angle methods were applied to the first vertical derivative of gravity data instead of Bouguer gravity data. We obtained the maps that show the previous lineaments which were designated by considering the maximum contours of horizontal gradient, analytic signal maps, and zero values of tilt angle, hyperbolic of tilt angle maps. The geometry of basement interface was also modeled utilizing the Parker–Oldenburg algorithm to understand the sediment thickness and coherency or incoherency between the gravity values and basement topography. The lineaments were held a candle to most current tectonic structure map of the study area. It was seen that the techniques used in this study are very effective to determine the old and new lineaments in the Shamakhy–Gobustan and Absheron regions. The epicenter distribution of earthquakes within the study area supports the new lineaments which are extracted by our interpretation. We concluded that better comprehension of Azerbaijan geostructures and its effect on the large scale works will be provided by means of this study.     

Forward plane-wave electromagnetic model in three dimensions using hybrid finite volume–integral equation scheme

We present a concept of the hybrid finite volume–integral equation technique for solving Maxwell's equation in a quasi-static form. The divergence correction was incorporated to improve the convergence and stability of the governing linear system equations which pose a challenge on the discretization of the curl–curl Helmholtz equation. A staggered finite volume approach is applied for discretizing the system of equations on a structured mesh and solved in a secondary field technique. The bi-conjugate gradient stabilizer was utilized with block incomplete lower-upper factorization preconditioner to solve the system of equation. To obtain the electric and magnetic fields at the receivers, we use the integral Green tensor scheme. We verify the strength of our hybrid technique with benchmark models relative to other numerical algorithms. Importantly, from the tested models, our scheme was in close agreement with the semi-analytical solution. It also revealed that the use of a quasi-analytical boundary condition helps to minimize the runtime for the linear system equation. Furthermore, the integral Green tensor approach to compute at the receivers demonstrates better accuracy compared with the conventional interpolation method. This adopted technique can be applied efficiently to the inversion procedure.     

Zircon U–Pb dating using LA-ICP-MS: Quaternary tephras in Boso Peninsula,Japan

Zircon U–Pb dating using LA-ICP-MS was applied to six Quaternary tephras in Boso Peninsula, central Japan: J1, Ks4, Ks5, Ks10, Ks11, and Ch2 in descending order. Accurate age determination of these tephras is of critical importance because they are widespread tephras in Japan and also relevant to a candidate site for the global boundary stratotype section and point of the early–middle Pleistocene boundary. Twenty grains were dated for each tephra and the following results were obtained. The J1 tephra had only 5 grains that yielded <2 Ma. The obtained age was ∼0.2 m.y. older than the stratigraphic age. No Quaternary ages were obtained from the Ks4 tephra. The Ks5 and Ks10 tephras had 10–12 grains that were ∼0.1–0.3 m.y. older than the stratigraphic age. The Ks11 tephra had 14 grains that yielded a weighted mean age of 0.52 ± 0.04 Ma (error reported as 95% confidence level), which was in agreement with the stratigraphic age. The Ch2 tephra had 16 grains that yielded a weighted mean age of 0.61 ± 0.02 Ma, which was also in agreement with the stratigraphic age. The good agreement between zircon U–Pb ages and the stratigraphy for Ks11 and Ch2 tephras validates the reliability of the established stratigraphy and our dating approach. The other tephras that yielded ∼0.1–0.3 m.y. older ages than the stratigraphy may indicate that the analyzed zircons were antecrysts that crystallized before eruption or they were detrital zircons incorporated during deposition.     

Interaction of Ca~(2+) and soil humic acid characterized by a joint experimental platform of potentiometric titration,UV–visible spectroscopy,and fluorescence spectroscopy

Rocky desertification has become a major environmental issue in the karst region of southwestern China.Karst rocky desertification was more severe in regions of limestone soil than in adjacent regions of other soils,despite the relatively higher soil organic matter(SOM)content in limestone soil. The underlying mechanism remains ambiguous. We speculated that the geochemical characteristics of limestone soils in the karst region plays an essential role, especially the high calcium content of limestone soil. To test this hypothesis, we collected limestone soil samples from a limestone soil profile in the southwestern China karst region and extracted humic acid(HA) from these limestone soil samples. We investigated the interaction of Ca~(2+)and three HA samples on a joint experimental platform, which consists of an automatic potentiometric titrator, a UV–visible spectrometer, and a Fluorescence spectrometer. HA solutions were titrated by Ca~(2+)and optical spectra of the HA solutions were monitored during the titration experiments. The results indicated that:(1) the interaction of Ca~(2+)and HA is a combined process of adsorption and complexation. Adsorption dominated the overall distribution behavior of Ca~(2+), which could be fit by Langmuir and Freundlich isotherm models.Complexation was distinguished only when the concentration of Ca~(2+)is low;(2) the changes of UV–visible spectroscopy and excitation–emission matrix fluorescencespectroscopy spectra of HA samples when they were binding with Ca~(2+)implied the apparent molecular size and structure of HA became larger and more complex;(3) the combination of Ca~(2+)and HA plays an important role in the SOM preservation of limestone soils but the stability of the Ca–HA association was relatively weak. The present study draws attention to maintaining the relatively higher Ca~(2+)concentration in limestone soils in ecologic restoration attempts in karst regions.     

Effects of temporally correlated infiltration on water flow in an unsaturated–saturated system

Effects of temporally correlated infiltration on water flow in an unsaturated–saturated system were investigated. Both white noise and exponentially correlated infiltration processes were considered. The moment equations of the pressure head (ψ) were solved numerically to obtain the variance and autocorrelation functions of ψ at 14 observation points. Monte Carlo simulations were conducted to verify the numerical results and to estimate the power spectrum of ψ (S _ψψ ). It was found that as the water flows through the system, the variance of the ψ (\( \sigma_{\psi }^{2} \)) were damped by the system: the deeper in the system, the smaller the \( \sigma_{\psi }^{2} \), and the larger the correlation timescale of the infiltration process (λ _I ), the larger the \( \sigma_{\psi }^{2} \). The unsaturated–saturated system gradually filters out the short-term fluctuations of ψ and the damping effect is most significant in the upper part of the system. The fluctuations of ψ is non-stationary at early time and becomes stationary as time progresses: the larger the value of λ _I , the longer the non-stationary period. The correlation timescale of the ψ (λ _ψ ) increases with depth and approaches a constant value at depth: the larger the value of λ _I , the larger the value of λ _ψ . The results of the estimated S _ψψ is consistent with those of the variance and autocorrelation function.