Bias and uncertainty in regression-calibrated models of groundwater flow in heterogeneous media

Groundwater models need to account for detailed but generally unknown spatial variability (heterogeneity) of the hydrogeologic model inputs. To address this problem we replace the large, m-dimensional stochastic vector β that reflects both small and large scales of heterogeneity in the inputs by a lumped or smoothed m-dimensional approximation γθ_∗, where γ is an interpolation matrix and θ_∗ is a stochastic vector of parameters. Vector θ_∗ has small enough dimension to allow its estimation with the available data. The consequence of the replacement is that model function f(γθ_∗) written in terms of the approximate inputs is in error with respect to the same model function written in terms of β, f(β), which is assumed to be nearly exact. The difference f(β) − f(γθ_∗), termed model error, is spatially correlated, generates prediction biases, and causes standard confidence and prediction intervals to be too small. Model error is accounted for in the weighted nonlinear regression methodology developed to estimate θ_∗ and assess model uncertainties by incorporating the second-moment matrix of the model errors into the weight matrix. Techniques developed by statisticians to analyze classical nonlinear regression methods are extended to analyze the revised method. The analysis develops analytical expressions for bias terms reflecting the interaction of model nonlinearity and model error, for correction factors needed to adjust the sizes of confidence and prediction intervals for this interaction, and for correction factors needed to adjust the sizes of confidence and prediction intervals for possible use of a diagonal weight matrix in place of the correct one. If terms expressing the degree of intrinsic nonlinearity for f(β) and f(γθ_∗) are small, then most of the biases are small and the correction factors are reduced in magnitude. Biases, correction factors, and confidence and prediction intervals were obtained for a test problem for which model error is large to test robustness of the methodology. Numerical results conform with the theoretical analysis.     

Multi-resolution adaptive modeling of groundwater flow and transport problems

Many groundwater flow and transport problems, especially those with sharp fronts, narrow transition zones, layers and fingers, require extensive computational resources. In this paper, we present a novel multi-resolution adaptive Fup approach to solve the above mentioned problems. Our numerical procedure is the Adaptive Fup Collocation Method (AFCM), based on Fup basis functions and designed through a method of lines (MOL). Fup basis functions are localized and infinitely differentiable functions with compact support and are related to more standard choices such as splines or wavelets. This method enables the adaptive multi-resolution approach to solve problems with different spatial and temporal scales with a desired level of accuracy using the entire family of Fup basis functions. In addition, the utilized collocation algorithm enables the mesh free approach with consistent velocity approximation and flux continuity due to properties of the Fup basis functions. The introduced numerical procedure was tested and verified by a few characteristic groundwater flow and transport problems, the Buckley–Leverett multiphase flow problem, the 1-D vertical density driven problem and the standard 2-D seawater intrusion benchmark–Henry problem. The results demonstrate that the method is robust and efficient particularly when describing sharp fronts and narrow transition zones changing in space and time.     

Stochastic modeling of Iranian earthquakes and estimation of ground motion for future earthquakes in Greater Tehran

Strong-motion data from eight significant well-documented earthquakes in Iran have been simulated using a stochastic modeling technique for finite faults proposed by Beresnev and Atkinson [Bull Seismol Soc Am 87 (1997) 67–84; Seism Res Lett 69 (1998) 27–32]. The database consists of 61 three-component records from eight earthquakes of magnitude ranging from M 6.3 to M 7.4, recorded at hypocentral distances up to 200 km. The model predictions are in good agreement with available Iranian strong-motion data as evidenced by near-zero average of differences between logarithms of the observed and predicted values for all frequencies. The strength factor, sfact, a quantity that controls the high-frequency radiation from the source is determined, on an event-by-event basis, by fitting simulated to observed response spectra.     

Upscaling unsaturated hydraulic parameters for flow through heterogeneous anisotropic sediments

We compare two methods for determining the upscaled water characteristics and saturation-dependent anisotropy in unsaturated hydraulic conductivity from a field-scale injection test. In both approaches an effective medium approximation is used to reduce a porous medium of M textures to an equivalent homogenous medium. The first approach is a phenomenological approach based on homogenization and assumes that moisture-based Richards’ equation can be treated like the convective–dispersive equation (CDE). The gravity term, dK_z(θ)/d(θ), analogous to the vertical convective velocity in the CDE, is determined from the temporal evolution of the plume centroid along the vertical coordinate allowing calculation of an upscaled K_z(θ). As with the dispersion tensor in the CDE, the rate of change of the second spatial moment in 3D space is used to calculate the water diffusivity tensor, D(θ), from which an upscaled K(θ) is calculated. The second approach uses the combined parameter scale inverse technique (CPSIT). Parameter scaling is used first to reduce the number of parameters to be estimated by a factor M. Upscaled parameters are then optimized by inverse modeling to produce an upscaled K(θ) characterized by a pore tortuosity–connectivity tensor, L. Parameters for individual textures are finally determined from the optimized parameters by inverse scaling using scale factors determined a priori. Both methods produced upscaled K(θ) that showed evidence of saturation dependent anisotropy. Flow predictions with the STOMP simulator, parameterized with upscaled parameters, were compared with field observations. Predictions based on the homogenization method were able to capture the mean plume behavior but could not reproduce the asymmetry caused by heterogeneity and lateral spreading. The CPSIT method captured the effects of heterogeneity and anisotropy and reduced the mean squared residual by nearly 90% compared to local-scale and upscaled parameters from the homogenization method. The Pacific Northwest National Laboratory is operated for the US Department of Energy by Battelle under Contract DE-AC05-76RL01830.     

Effect of zonal E × B plasma drift on electron density in the low-latitude ionospheric <Emphasis Type="Italic">F</Emphasis> region at a solar activity maximum near vernal equinox

The variations in the density of the ionospheric F2 layer maximum (NmF2) under the action of the zonal plasma drift perpendicularly to the magnetic (B) and electric (E) fields in the direction geomagnetic west-geomagnetic east have been studied using the three-dimensional nonstationary theoretical model of electron and ion densities (N _e and N _i ) and temperatures (T _e and T _i ) in the low-latitude and midlatitude ionospheric F region and plasmasphere. The method of numerical calculations of N _e , N _i , T _e , and T _i , including the advantages of the Lagrangian and Eulerian methods, is used in the model. A dipole approximation of the geomagnetic field (B), taking into account the non-coincidence of the geographic and geomagnetic poles and differences between the positions of the Earth’s and geomagnetic dipole centers, is accepted in the calculations. The calculated NmF2 and altitudes of the F2 layer maximum (hmF2) have been compared with these quantities measured at 16 low-latitude ionospheric sounding stations during the geomagnetically quiet period October 11–12, 1958. This comparison made it possible to correct the input model parameters: the NRLMSISE-00 model [O], the meridional component of the neutral wind velocity according to the HWW90 model, and the meridional component of the equatorial plasma drift due to the electric field specified by the empirical model. It has been indicated that the effect of the zonal E × B plasma drift on NmF2 can be neglected under daytime conditions and changes in NmF2 and hmF2 under the action of this drift are insignificant under nighttime conditions north of 25° and south of ?26° geomagnetic latitude. The effect of the zonal E × B plasma drift on NmF2 and hmF2 is most substantial in the nightside ionosphere approximately from ?20° to 20° geomagnetic latitude, and the neglect of this drift results in an up to 2.4-fold underestimation of NmF2. The found dependence of the effect of the zonal E × B plasma drift on NmF2 and hmF2 on geomagnetic latitude is related to the longitudinal asymmetry of B, asymmetry of the neutral wind about the geomagnetic equator, and changes in the meridional E × B plasma drift at a change in geomagnetic longitude.     

Refraction of the principal stress trajectories and the stress jumps on faults and contact surfaces: Part 1. Non-constrained regular trajectories

Rock masses contain ubiquitous multiscale heterogeneities, which (or whose boundaries) serve as the surfaces of discontinuity for some characteristics of the stress state, e.g., for the orientation of principal stress axes. Revealing the regularities that control these discontinuities is a key to understanding the processes taking place at the boundaries of the heterogeneities and for designing the correct procedures for reconstructing and theoretical modeling of tectonic stresses. In the present study, the local laws describing the refraction of the axes of extreme principal stresses T ₁ (maximal tension in the deviatoric sense) and T ₃ (maximal compression) of the Cauchy stress tensor at the transition over the elementary area n of discontinuity whose orientation is specified by the unit normal n are derived. It is assumed that on the area n of discontinuity, frictional contact takes place. No hypotheses are made on the constitutive equations, and a priori constraints are not posed on the orientation on the stress axes. Two domains, which adjoin area n on the opposite sides and are conventionally marked + and ?, are distinguished. In the case of the two-dimensional (2D) stress state, any principal stress axis on passing from domain ? to domain + remains in the same quadrant of the plane as the continuation of this axis in domain +. The sign and size of the refraction angle depend on the sign and amplitude of the jump of the normal stress, which is tangential to the surface of discontinuity. In the three-dimensional (3D) case, the refraction of axes T ₁ and T ₃ should be analyzed simultaneously. For each side, + and ?, the projections of the T ₁ and T ₃ axes on the generally oriented plane n form the shear sectors S ⁺ and S ^?, which are determined unambiguously and to whose angular domains the possible directions p ⁺ and p ^? of the shear stress vectors belong. In order for the extreme stress axes T ₁ ⁺ ,T ₃ ⁺ and T ₁ ^? , T ₃ ^? to be statically compatible on the generally oriented plane n, it is required that sectors S ⁺ and S ^? had a nonempty intersection. The direction vectors p ⁺ and p ^? are determined uniquely if, besides axes T ₁ ^? , T ₃ ^? and T ₁ ⁺ , T ₃ ⁺ , also the ratios of differential stresses R ⁺ and R ^? (0 ≤ R ^± ≤ 1) are known. This is equivalent to specifying the reduced stress tensors T _R ⁺ and T _R ^? The necessary condition for tensors T _R ⁺ and T _R ^? being statically compatible on plane n is the equality p ⁺ = p ^?. In this paper, simple methods are suggested for solving the inverse problem of constructing the set of the orientations of the extreme stress axes from the known direction p of the shear stress vector on plane n and from the data on the shear sector. Based on these methods and using the necessary conditions of local equilibrium on plane n formulated above, all the possible orientations of axes T ₁ ⁺ , T ₃ ⁺ are determined if the projections of axes T ₁ ^? , T ₃ ^? axes on side — are given. The angle between the projections of axes T ₁ ⁺ , T ₁ ^? and/or T ₃ ⁺ , T ₃ ^? on the plane can attain 90°. Besides the general case, also the particular cases of the contact between the degenerate stress states and the special position of plane n relative to the principal stress axes are thoroughly examined. Generalization of the obtained results makes it possible to plot the local diagram of the orientations of axes T ₁ ⁺ , T ₃ ⁺ for a given sector S ^?. This diagram is a so-called stress orientation sphere, which is subdivided into three pairs of areas (compression, tension, and compression-extension). The tension and compression zones cannot contain the poles of T ₃ ⁺ and T ₁ ⁺ axes, respectively. The compression-extension zones can contain the poles of either T ₁ ⁺ or T ₃ ⁺ axis but not both poles simultaneously. In the particular case when the shear stress vector has a unique direction p ^? on side ?, the areas of compression-extension disappear and the diagram is reduced to a beach-ball plot, which visualizes the focal mechanism solution of an earthquake. If area n is a generally oriented plane and if the orientation of the pairs of the statically compatible axes T ₁ ^? , T ₃ ^? and T ₁ ⁺ , T ₃ ⁺ is specified, then, the stress values on side + are uniquely determined from the known stress values on side ?. From the value of differential stress ratio R ^?, one can calculate the value of R ⁺, and using the values of the principal stresses on side ?, determine the total stress tensor T ⁺ on side +. The obtained results are supported by the laboratory experiments and drilling data. In particular, these results disclose the drawbacks of some established notions and methods in which the possible refraction of the stress axes is unreasonably ignored or taken into account improperly. For example, it is generally misleading to associate the slip on the preexisting fault with the orientation of any particular trihedron of the principal stress axes. The reconstruction should address the potentially statically compatible principal stress axes, which are differently oriented on opposite sides of the fault plane. The fact that, based on the orientation of the intraplate principal stresses at the base of the lithosphere, one cannot make a conclusion on the active or passive influence of the mantle flows on the lithospheric plate motion is another example. The present relationships linking the stress values on the opposite sides of the fault plane on which the orientations of the principal stress axes are known demonstrate the incorrectness of the existing methods, in which the reduced stress tensors within the material domains are reconstructed without allowance for the dynamic interaction of these domains with their neighbors. In addition, using the obtained results, one can generalize the notion of the zone of dynamical control of a fault onto the case of the existence of discontinuities in this region and analyze the stress transfer across the system of the faults.     

Mechanism of formation of <Emphasis Type="Italic">Q</Emphasis>-disturbances in the <Emphasis Type="Italic">F</Emphasis><Subscript>2</Subscript> region of the equatorial ionosphere

Using model simulations, the morphological picture (revealed earlier) of the disturbances in the F ₂ region of the equatorial ionosphere under quiet geomagnetic conditions (Q-disturbances) is interpreted. It is shown that the observed variations in the velocity of the vertical E × B plasma drift, related to the zonal E _y component of the electric field, are responsible for the formation of Q-disturbances. The plasma recombination at altitudes of the lower part of the F ₂ region and the dependence of the rate of this process on heliogeophysical conditions compose the mechanism of Q-disturbance formation at night. The daytime positive Q-disturbances are caused exclusively by a decrease in the upward E × B drift, and this type of disturbances could be related to the known phenomenon of counter electrojet. Possible causes of formation of the daytime negative Q-disturbances are discussed.     

An efficient,high-order probabilistic collocation method on sparse grids for three-dimensional flow and solute transport in randomly heterogeneous porous media

In this study, a probabilistic collocation method (PCM) on sparse grids is used to solve stochastic equations describing flow and transport in three-dimensional, saturated, randomly heterogeneous porous media. The Karhunen–Loève decomposition is used to represent log hydraulic conductivity Y=ln_Ks$Y=\ln K_s$. The hydraulic head h   and average pore-velocity v$\mathbf{v}$ are obtained by solving the continuity equation coupled with Darcy’s law with random hydraulic conductivity field. The concentration is computed by solving a stochastic advection–dispersion equation with stochastic average pore-velocity v$\mathbf{v}$ computed from Darcy’s law. The PCM approach is an extension of the generalized polynomial chaos (gPC) that couples gPC with probabilistic collocation. By using sparse grid points in sample space rather than standard grids based on full tensor products, the PCM approach becomes much more efficient when applied to random processes with a large number of random dimensions. Monte Carlo (MC) simulations have also been conducted to verify accuracy of the PCM approach and to demonstrate that the PCM approach is computationally more efficient than MC simulations. The numerical examples demonstrate that the PCM approach on sparse grids can efficiently simulate solute transport in randomly heterogeneous porous media with large variances.     

Efficient flow and transport simulations in reconstructed 3D pore geometries

Upscaling pore-scale processes into macroscopic quantities such as hydrodynamic dispersion is still not a straightforward matter for porous media with complex pore space geometries. Recently it has become possible to obtain very realistic 3D geometries for the pore system of real rocks using either numerical reconstruction or micro-CT measurements. In this work, we present a finite element–finite volume simulation method for modeling single-phase fluid flow and solute transport in experimentally obtained 3D pore geometries. Algebraic multigrid techniques and parallelization allow us to solve the Stokes and advection–diffusion equations on large meshes with several millions of elements. We apply this method in a proof-of-concept study of a digitized Fontainebleau sandstone sample. We use the calculated velocity to simulate pore-scale solute transport and diffusion. From this, we are able to calculate the a priori emergent macroscopic hydrodynamic dispersion coefficient of the porous medium for a given molecular diffusion D_m of the solute species. By performing this calculation at a range of flow rates, we can correctly predict all of the observed flow regimes from diffusion dominated to convection dominated.     

Inverse modeling of dissolved O2 and δ18O-DO to estimate aquatic metabolism,reaeration and respiration isotopic fractionation: effects of variable light regimes and input uncertainties

A transient model, hereafter referred to as ROM-TM, was developed to quantify river ecosystem metabolic rates and reaeration rates from field observation of changes in dissolved O₂ (DO) and the ratio of ¹⁸O to ¹⁶O in DO (δ¹⁸O-DO). ROM-TM applies an inverse modeling approach and is programmed using MATLAB. Parameters describing photosynthesis, ecosystem respiration, gas exchange, and isotopic fractionation, such as maximum photosynthetic rate (P _m ), photosynthetic efficiency parameter (a), respiration rate at 20 °C (R ₂₀ ), gas exchange coefficient (K), respiration isotopic fractionation factor (a _R ), and photorespiration coefficient (β _R ), can be abstracted by minimizing the sum of square errors between the fitted data and the observed field data. Then DO and δ¹⁸O-DO time series can be reconstructed using estimated parameters and input variables. Besides being capable of teasing apart metabolic processes and gas exchange to provide daily average estimates of metabolic parameters at the ecosystem scale, ROM-TM can be used to address issues related to light including light saturation phenomena at the ecosystem level, the effect of cloud cover on the metabolic balance, and photorespiration. Error and uncertainty analysis demonstrates that ROM-TM is stable and robust for the random errors of DO time series. The photosynthetic parameters P _m and a are more sensitive than other parameters to lower-resolution time series data.     

A coupled model for simulating surface water and groundwater interactions in coastal wetlands

Coastal wetlands are characterized by strong, dynamic interactions between surface water and groundwater. This paper presents a coupled model that simulates interacting surface water and groundwater flow and solute transport processes in these wetlands. The coupled model is based on two existing (sub) models for surface water and groundwater, respectively: ELCIRC (a three‐dimensional (3‐D) finite‐volume/finite‐difference model for simulating shallow water flow and solute transport in rivers, estuaries and coastal seas) and SUTRA (a 3‐D finite‐element/finite‐difference model for simulating variably saturated, variable‐density fluid flow and solute transport in porous media). Both submodels, using compatible unstructured meshes, are coupled spatially at the common interface between the surface water and groundwater bodies. The surface water level and solute concentrations computed by the ELCIRC model are used to determine the boundary conditions of the SUTRA‐based groundwater model at the interface. In turn, the groundwater model provides water and solute fluxes as inputs for the continuity equations of surface water flow and solute transport to account for the mass exchange across the interface. Additionally, flux from the seepage face was routed instantaneously to the nearest surface water cell according to the local sediment surface slope. With an external coupling approach, these two submodels run in parallel using time steps of different sizes. The time step (Δt_g) for the groundwater model is set to be larger than that (Δt_s) used by the surface water model for computational efficiency: Δt_g = M × Δt_s where M is an integer greater than 1. Data exchange takes place between the two submodels through a common database at synchronized times (e.g. end of each Δt_g). The coupled model was validated against two previously reported experiments on surface water and groundwater interactions in coastal lagoons. The results suggest that the model represents well the interacting surface water and groundwater flow and solute transport processes in the lagoons. Copyright © 2011 John Wiley & Sons, Ltd.     

Numerical experiments of air–ice drag coefficient and its impact on ice–ocean coupled system in the Sea of Okhotsk

Sensitivity study of the air–ice drag coefficient C _Dai is presented with an ice–ocean coupled model for the Sea of Okhotsk. The C _Dai?×?10³ value is varied from 2 to 5 based on the direct measurements in the region. The maximum volume transport of the East Sakhalin Current and the mean sea ice velocity were intensified as C _Dai increased. The sensitivity experiment with the ice–water drag coefficient C _Diw showed that the East Sakhalin Current volume transport is hardly affected by C _Diw but significantly intensified by C _Dai. While the ice drift in the off-ice-edge direction was intensified by the increase in C _Dai and the decrease in C _Diw, the ice edge location was nearly unchanged. This was due to melting caused by the relatively warm water inflow from the North Pacific. That is, sea ice extent in the region is strongly influenced by melting caused by a large ice–ocean heat transfer. In the active melting regions, the ice–ocean heat transfer of more than 100 W/m² occurred even in mid-winter. This is the same order as the cooling by air in winter, and a heat insulation capacity of sea ice is weakened in such regions.     

Geomagnetically induced currents in power lines according to data on geomagnetic variations

The time derivative (d H/dt) of the geomagnetic field horizontal component (H) for seven intervals of high geomagnetic activity in 2003–2005 has been calculated, based on the data of Alma-Ata, Novosibirsk, and Irkutsk observatories, in order to estimate the probability of appearance of geomagnetically induced current (GIC), the value of which is linearly dependent on d H/dt, in power lines on the territory of Kazakhstan. The distributions of the H and d H/dt directions have been constructed; in this case it was most interesting that these distributions were narrow and extended along the magnetic meridian for Alma-Ata and were wider angular for Novosibirsk and Irkutsk. It has been indicated that large H values, determining significant GIC values, took place at a sudden commencement of strong storms, which had a character of a pulsed disturbance of the geomagnetic field, and during large-amplitude geomagnetic field pulsations. The duration of the periods of large |d H/dt| values exceeding the threshold (30 nT/min), when GIC could cause unwanted consequences in power circuits, has been determined.     

Correlation between abnormal trends in the spontaneous fields of tectonic plates and strong seismicities

Tectonic activities, electrical structures, and electromagnetic environments are major factors that affect the stability of spontaneous fields. The method of correlating regional synchronization contrasts(CRSC) can determine the reliability of multi-site data trends or shortimpending anomalies. From 2008 to 2013, there were three strong earthquake cluster periods in the North–South seismic belt that lasted for 8–12 months. By applying the CRSC method to analyze the spontaneous field E_(SP) at 25 sites of the region in the past 6 years, it was discovered that for each strong earthquake cluster period, the E_(SP) strength of credible anomalous trends was present at minimum 30%of the stations. In the southern section of the Tan-Lu fault zone, the E_(SP) at four main geoelectric field stations showed significant anomalous trends after June 2015, which could be associated with the major earthquakes of the East China Sea waters(MS7.2) in November 2015 and Japan's Kyushu island(MS7.3) in April 2016.     

The relation between solar activity and orientation of the solar wind electric field relative to the Earth’s magnetic moment

The relation of the Kp index of geomagnetic activity to the solar wind electric field (E _SW) and the projection of this field onto the geomagnetic dipole has been estimated. An analysis indicated that the southward component of the IMF vector (B _z < 0) is the main geoeffective parameter, as was repeatedly indicated by many researchers. The presence of this component in any combinations of the interplanetary medium parameters is responsible for a high correlation between such combinations and geomagnetic activity referred to by the authors of different studies. Precisely this field component also plays the main role in the relation between the Kp index and the relative orientation of E _SW and the Earth’ magnetic moment.