Real-time multipurpose reservoir operation: a case study

Abstract

A real-time operational methodology has been developed for multipurpose reservoir operation for irrigation and hydropower generation with application to the Bhadra reservoir system in the state of Karnataka, India. The methodology consists of three phases of computer modelling. In the first phase, the optimal release policy for a given initial storage and inflow is determined using a stochastic dynamic programming (SDP) model. Streamflow forecasting using an adaptive AutoRegressive Integrated Moving Average (ARIMA) model constitutes the second phase. A real-time simulation model is developed in the third phase using the forecast inflows of phase 2 and the operating policy of phase 1. A comparison of the optimal monthly real-time operation with the historical operation demonstrates the relevance, applicability and the relative advantage of the proposed methodology.     

Performance evaluation of implicit stochastic reservoir operation optimization supported by long-term mean inflow forecast

This study applies implicit stochastic optimization (ISO) to develop monthly operating rules for a reservoir located in Northeast Brazil. The proposed model differs from typical ISO applications as it uses the forecast of the mean inflow for a future horizon instead of the current-month inflow. Initially, a hundred different 100-year monthly inflow scenarios are synthetically generated and employed as input to a deterministic operation optimization model in order to build a database of optimal operating data. Later, such database is used to fit monthly reservoir rule curves by means of nonlinear regression analysis. Finally, the established rule curves are validated by operating the system under 100 new inflow ensembles. The performance of the proposed technique is compared with those provided by the standard reservoir operating policy (SOP), stochastic dynamic programming (SDP) and perfect-forecast deterministic optimization (PFDO). Different forecasting horizons are tested. For all of them, the results indicate the feasibility of using ISO in view of its lower vulnerability in contrast to the SOP as well as the proximity of its operations with those by PFDO. The results also reveal that there is an optimal choice for the forecasting horizon. The comparison between ISO and SDP shows small differences between both, justifying the adoption of ISO for its simplified mathematics as opposed to SDP.     

A stability analysis of a conduit flow model for lava dome eruptions

Periodic variations in magma discharge rate and ground deformation have been commonly observed during lava dome eruptions. We performed a stability analysis of a conduit flow model by Barmin et al. [Barmin, A., Melnik, O., Sparks, R.S.J., 2002. Periodic behavior in lava dome eruptions. Earth and Planetary Science Letters 199 (1-2), 173–184], in which the periodic variations in magma flow rate and chamber pressure are reproduced as a result of the temporal and spatial changes of the magma viscosity controlled by the kinetics of crystallization. The model is reduced to a dynamical system where the time derivatives of the magma flow rate (dQ/dt) and the chamber pressure (dP/dt) are functions of Q and P evaluated at a shifted time t − t_?. Here, the time delay t_? represents the time for the viscosity of fluid particle to increase in a conduit. The dynamical system with time delay is approximated by a simple two-dimensional dynamical system of Q and P where t_? is given as a parameter. The results of our linear stability analyses for these dynamical systems indicate that the transition from steady to periodic flow depends on nonlinearities in the steady state relation between Q and P. The steady state relation shows a sigmoidal curve in Q − P phase plane; its slope has negative values at intermediate flow rates. The steady state solutions become unstable, and hence P and Q oscillate periodically, when the negative slope of the steady state relation ([dP/dQ]_S) exceeds a critical value; that is [dP/dQ]_S < − t_?γ/(2V_ch), where V_ch is the chamber volume and γ is an elastic constant which is related to the rigidity of chamber wall. We also found that the period and the pattern of oscillation of the conduit flow primarily depend on a quantity defined by LV_ch/r⁴, where L is the conduit length and r is the conduit radius.     

Hydroclimatology of the Volta River Basin in West Africa: Trends and variability from 1901 to 2002

Long-term historical records of rainfall (P), runoff (Q) and other climatic factors were used to investigate hydrological variability and trends in the Volta River Basin over the period 1901-2002. Potential (E_p) and actual evaporation (E), rainfall variability index (δ), Budyko’s aridity index (I_A), evaporation ratio (C_E) and runoff ratio (C_Q) were estimated from the available hydroclimatological records. Mann-Kendall trend analysis and non-parametric Sen’s slope estimates were performed on the respective time series variables to detect monotonic trend direction and magnitude of change over time.Rainfall variability index showed that 1968 was the wettest year (δ = +1.75) while 1983 was the driest (δ = −3.03), with the last three decades being drier than any other comparable period in the hydrological history of the Volta. An increase of 0.2 mm/yr² (P < 0.05) was observed in E_p for the 1901-1969 sub-series while an increased of 1.8 mm/yr² (P < 0.01) was recorded since 1970. Rainfall increased at the rate of 0.7 mm/yr² or 49 mm/yr between 1901 and 1969, whereas a decrease of 0.2 mm/yr² (6 mm/yr) was estimated for 1970-2002 sub-series. Runoff increased significantly at the rate of 0.8 mm/yr (23 mm/yr) since 1970. Runoff before dam construction was higher (87.5 mm/yr) and more varied (CV = 41.5%) than the post-dam period with value of 73.5 mm/yr (CV = 23.9%). A 10% relative decrease in P resulted in a 16% decrease in Q between 1936 and 1998. Since 1970, all the months showed increasing runoff trends with significant slopes (P < 0.05) in 9 out of the 12 months. Possible causes, such as climate change and land cover change, on the detected changes in hydroclimatology are briefly discussed.     

Yield and overflow from a finite stochastic reservoir

The paper discusses the overflow (spillage) and yield rates and the total overflow and total yield over a specified time from a finite discrete stochastic reservoir, in which the yieldY _t during the working interval (t,t+1) is a function of the storageZ _t at timet, the inflow sequence {X _t } being IID.The distribution vector of the spillage rate at timet is a telescoped version of the distribution of a certain Markovian variable whose transition matrix is derived. Formulae are given for the distribution of the total spillageW _h given suitable initial conditions, forh=1,2,3; and a simple expression derived forE(W _h ).     

Information,artifacts, and noise in dQ/dt − Q recession analysis

The plotting of the time rate of change in discharge dQ/dt versus discharge Q has become a widely used tool for analyzing recession data since Brutseart and Nieber [Water Resour Res 13 (1977) 637–643] proposed the method. Typically the time increment Δt over which the recession slope dQ/dt is approximated is held constant. It is shown here this that leads to upper and lower envelopes in graphs of log(−dQ/dt) versus log(Q) that have been observed in previous studies but are artifacts. The use of constant time increments also limits accurate representation of the recession relationship to the portion of the hydrograph for which the chosen time increment is appropriate. Where dQ/dt varies by orders of magnitude during recession, this may exclude much of the hydrograph from analysis. In response, a new method is proposed in which Δt for each observation in time is properly scaled to the observed drop in discharge ΔQ. It is shown, with examples, how the new method can succeed in exposing the underlying relationship between dQ/dt and Q where the standard method fails.     

Coda-Q and its lapse time dependence analysis in the interaction zone of the Dinarides,the Alps and the Pannonian basin

The single backscattering model was used to estimate total attenuation of coda waves (Q_c) of local earthquakes recorded on eight seismological stations in the complex area of the western continental Croatia. We estimated Q₀ and n, parameters of the frequency dependent coda-Q using the relation Q_c = Q₀fⁿ. Lapse time dependence of these parameters was studied using a constant 30 s long time window that was slid along the coda of seismograms. Obtained Q_c were distributed into classes according to their lapse time, t_L. For t_L = 20–50 s we estimated Q₀ = 45–184 and n = 0.49–0.94, and for t_L = 60–100 s we obtained Q₀ = 119–316 and n = 0.37–0.82. There is a tendency of decrease of parameter n with increasing Q₀, and vice versa. The rates of change of both Q₀ and n seem to decrease for lapse times larger than 50–80 s, indicating an alteration in rock properties controlling coda attenuation at depths of about 100–160 km. A very good correlation was found between the frequency dependence parameter n and the Moho depths for lapse times of 50, 60 and 70 s.     

Scattering and intrinsic attenuation in Cairo metropolitan area using genetic algorithm

A total number of 46 local earthquakes (2.0≤ML≤4.0) recorded in the period 2000–2011 by the Egyptian seismographic network (ENSN) were used to estimate the total (Q_t⁻¹), intrinsic (Q_i⁻¹) and scattering attenuation (Q_sc⁻¹) in Cairo metropolitan area, Egypt. The multiple lapse time window analysis (MLTWA) under the assumption of multiple isotropic scattering with uniform distribution of scatters was firstly applied to estimate the pair of L_e⁻¹, the extinction length inverse, and B₀, the seismic albedo, in the frequency range 3–24 Hz. To take into account the effect of a depth-dependent earth model, the obtained values of B₀ and L_e⁻¹ were corrected for an earth structure characterized by a transparent upper mantle and a heterogeneous crust. The estimated values of Q_t⁻¹, Q_sc⁻¹ and Q_i⁻¹ exhibited frequency dependences. The average frequency-dependent relationships of attenuation characteristics estimated for the region are found to be: Q_t⁻¹=(0.015±0.008)f ^{(−1.02±0.02)}, Q_sc⁻¹=(0.006±0.001)f ^{(−1.01±0.02)}, and Q_i⁻¹=(0.009±0.008)f ^{(−1.03±0.02)}; showing a predominance of intrinsic absorption over scattering attenuation. This finding implies that the pore-fluid contents may have great effect on the attenuation mechanism in the upper crust where the River Nile is passing through the study area. The obtained results are comparable with those obtained in other tectonic regions.     

基于生物水文响应模型的刺激四大家鱼繁殖生态调度目标量化方法研究——以三峡水库为例

为了改善生态调度效果,提高生态调度实践成功率,采用医学领域适用性广、准确性高、直观性强的受试者工作特征曲线法(receiver operating characteristic, ROC),以三峡水库为例,提出了一套基于生物水文响应模型的刺激四大家鱼繁殖生态调度目标量化方法,包括涨水事件界定及生态水文指标选择、鱼类自然繁殖事件及有效繁殖响应界定、生物-水文响应模型构建与多指标判别优选等多个环节。运用三峡水库运行后2013—2019年(除2016年)6年47场次宜昌江段的涨水事件和24次四大家鱼有效繁殖响应的实际观测数据构建模型并量化三峡水库生态调度目标。结果表明:(1)依据初始流量Q₀、洪峰流量Q_peak、流量总增长量Q_sum和流量日增长率dQ这4个指标所建立的单指标生物-水文响应模型具有较好判别四大家鱼自然繁殖对水文条件是否响应的能力;(2)各评价指标Q₀、Q_peak、Q_sum、dQ的最佳阈值分别为14960 m³/s、19610 m³/s、4050 m³/s及1405 m³/(s·d);单指标模型性能排序为dQ＞Q_peak＞Q_sum＞Q₀;(3)依据dQ或Q_peak的指标组合方式建立的多指标生物-水文响应模型综合性能最佳。结果表明,ROC模型方法结果形式简单、预见性强,据此量化生态调度目标,可为其他流域、水库的鱼类繁殖响应判别和水库生态调度方案制定提供参考借鉴。     

Attenuation characteristics of coda waves and estimation ofQc values in eastern China

For short-period near-earthquake records in eastern China, from the empirical attenuation formula of coda ground motion amplitudeA with timeτ: lgA=G?2. 235 lgτ, using the single scattering theory modified with epicentral distance, we obtain the curve family of corrected coda amplitudeA _c(r,t), andω/2Q _c values for each time interval of coda. From this,Q _c(f,h) values, which correspond to each observational average frequency and sampling depth, are calculated. The results substantially agree with those observationalQ _c values in Yunnan, Beijing and central Asia.     

Advance first order second moment (AFOSM) method for single reservoir operation reliability analysis: a case study

Reservoir system reliability is the ability of reservoir to perform its required functions under stated conditions for a specified period of time. In classical method of reservoir system reliability analysis, the operation policy is used in a simple simulation model, considering the historical/synthetic inflow series and a number of physical bounds on a reservoir system. This type of reliability analysis assumes a reservoir system as fully failed or functioning, called binary state assumption. A number of researchers from various research backgrounds have shown that the binary state assumption in the traditional reliability theory is not extensively acceptable. Our approach to tackle the present problem space is to implement the algorithm of advance first order second moment (AFOSM) method. In this new method, the inflow and reservoir storage are considered as uncertain variables. The mean, variance and covariance of uncertain variables are determined using moment values of reservoir state variables. For this purpose, a stochastic optimization model developed based on the constraint state formulation is applied. The proposed model of reliability analysis is used to a real case study in Iran. As a result, monthly probabilities of water allocation were computed from AFOSM method, and the outputs were compared with those from Monte Carlo method. The comparison shows that the outputs from AFOSM method are similar to those from the Monte Carlo method. In term of practical use of this study, the proposed method is appropriate to determine the monthly probability of failure in water allocation without the aid of simulation.     

Yield and overflow from a finite stochastic reservoir

The paper discusses the overflow (spillage) and yield rates and the total overflow and total yield over a specified time from a finite discrete stochastic reservoir, in which the yieldY _t during the working interval (t,t+1) is a function of the storageZ _t at timet, the inflow sequence {X _t } being IID.The distribution vector of the spillage rate at timet is a telescoped version of the distribution of a certain Markovian variable whose transition matrix is derived. Formulae are given for the distribution of the total spillageW _h given suitable initial conditions, forh=1,2,3; and a simple expression derived forE(W _h ).The distribution of the yield rateY _t is a trivial modification of the storage distribution. As for the total yieldR _t =Y ₁+...+Y _t , it is shown that the bivariate sequence {R _t ,Z _t } is first-order Markovian, whereZ _t is the storage at timet. The transition matrix of this process is obtained and the method of evaluating the marginal distribution of the total yieldR _t is exemplified.     

Exploring the nature of non-Fickian transport in laboratory experiments

Observations of non-Fickian transport in sandbox experiments [Levy M, Berkowitz B. Measurement and analysis of non-Fickian dispersion in heterogeneous porous media. J Contam Hydrol 2003;64:203–26] were analyzed previously using a power law tail ψ(t) ∼ t^−1−β with 0 < β < 2 for the spectrum of transition times comprising a tracer plume migration. For each sandbox medium a choice of β resulted in an excellent fit to the breakthrough curve (BTC) data, and the value of β decreased slowly with increasing flow velocity. Here, the data are reanalyzed with the full spectrum of ψ(t) gleaned from analytical calculations [Cortis A, Chen Y, Scher H, Berkowitz B. Quantitative characterization of pore-scale disorder effects on transport in “homogeneous” granular media. Phys Rev E 2004;10(70):041108. doi: 10.1103/PhysRevE.70.041108], numerical simulations [Bijeljic B, Blunt MJ. Pore-scale modeling and continuous time random walk analysis of dispersion in porous media. Water Resour Res 2006;42:W01202. doi: 10.1029/2005WR004578] and permeability fields [Di Donato G, Obi E-O, Blunt MJ. Anomalous transport in heterogeneous media demonstrated by streamline-based simulation. Geophys Res Lett 2003;30:1608–12s. doi: 10.1029/2003GL017196]. We represent the main features of the full spectrum of transition times with a truncated power law (TPL), ψ(t) ∼ (t₁ + t)^−1−βexp(−t/t₂), where t₁ and t₂ are the limits of the power law spectrum. An excellent fit to the entire BTC data set, including the changes in flow velocity, for each sandbox medium is obtained with a single set of values of t₁, β, t₂. The influence of the cutoff time t₂ is apparent even in the regime t < t₂. Significantly, we demonstrate that the previous apparent velocity dependence of β is a result of choosing a pure power law tail for ψ(t). The key is the change in the log–log slope of the TPL form of ψ(t) with a shifting observational time window caused by the change in the mean velocity. Hence, the use of the full spectrum of ψ(t) is not only necessary for the transition to Fickian behavior, but also to account for the dynamics of these laboratory observations of non-Fickian transport.     

Impact of forest degradation on streamflow regime and runoff response to rainfall in the Garhwal Himalaya,Northwest India

Baseflows have declined for decades in the Lesser Himalaya but the causes are still debated. This paper compares variations in streamflow response over three years for two similar headwater catchments in northwest India with largely undisturbed (Arnigad) and highly degraded (Bansigad) oak forest. Hydrograph analysis suggested no catchment leakage, thereby allowing meaningful comparisons. The mean annual runoff coefficient for Arnigad was 54% (range 44–61%) against 62% (53–69%) at Bansigad. Despite greater total runoff Q_t (by 250 mm year^–1), baseflow at Bansigad ceased by March, but was perennial at Arnigad (making up 90% of Q_t vs. 51% at Bansigad). Arnigad storm flows, Q_s, were modest (8–11% of Q_t) and occurred mostly during monsoons (78–98%), while Q_s at Bansigad was 49% of Q_t and occurred also during post-monsoon seasons. Our results underscore the importance of maintaining soil water retention capacity after forest removal to maintain baseflow levels.

EDITOR D. Koutsoyiannis; ASSOCIATE EDITOR D. Gerten     

Comparison of the attenuation properties for two different areas in eastern Anatolia,Turkey

In this study, the attenuation properties of the crust and the quality factor of S wave in eastern Anatolia (Turkey) were determined by local earthquakes for two different areas, Oltu and Erzurum. Seismic wave attenuation can be changed with high pressure or structural effects. Therefore, we argued that the estimation of attenuation coefficient in seismic active zones in Eastern Anatolia is a very useful tool to determine seismic activities. It uses regional waveform data set from two stations, OLT and ERZ, for 95 events that occurred in these regions between 2001 and 2005. The attenuation has been determined using the Chobra–Alexeev model based on the epicenter distance–amplitude relations. This model allows for investigation of the effects of variations in attenuation properties for different areas. We introduced a new magnitude formula for these areas using the amplitude normalization methods for reference values M_L=4, so as to correct effects of the magnitudes. We also determined velocity of seismic waves. The average attenuation coefficient (α), average quality factor (Q_s) and P and S waves velocities were obtained with normalized amplitude values for Erzurum (ERZ) and Oltu (OLT) as 0.0135 km⁻¹, 37, 6.20 km/s and 3.38 km/s and 0.0151, 34, 6.13 and 3.48.     

A scattering model for seismic attenuation and its global applications

The attenuation characteristics of Indian lithosphere and its comparison with different tectonic settings in the world are determined from the observations of the Q for L_g(Q_Lg)-, and S(Q_S)-waves in the 1-30 Hz frequency range. The scattering is approximated with a Gaussian distribution of spherical scatterers. To approximate single scattering, we use Dainty's [Geophy. Res. Lett. 8 (11) (1981) 1126] model that attenuation is given by 1/Q(ω) = 1/Q_i + g(ω)v/ω, where Q_i is intrinsic Q due to anelastic attenuation, v is shear wave velocity, ω is angular frequency, g = ∫n(a)σ da is the total scattering coefficient for S-to-S scattering, n(a) da is the number of scattering spheres of radius a per unit volume, and σ is the scattering cross-section for the sphere. We find that if n(a) is described by a simple two parameter (a₀ and c) Gaussian of amplitude c and standard deviation and mean a₀, the attenuation data for different regions of the world are well approximated over the frequency band of seismic observations. Our major findings are: (1) the maximum effect of scattering on attenuation occurs at 0.84 Hz or a wavelength of 4.16 km; (2) the values of g are frequency dependent. Values of g are of the order of 10⁻³ km⁻¹ at 1-30 Hz, varying from 0.0031 to 0.01 and 0.001 to 0.0083 km⁻¹ for tectonically active and stable regions, respectively; (3) regions of active tectonics and seismicity generally have lower Q_i values (1000) than that in stable regions (2000); and (4) regions of high Q_i value exhibit low intensity of scattering.     

Attenuation of Broadband P and S Waves in Tonga: Observations of Frequency Dependent Q

—Teleseismic broadband recordings of intermediate and deep focus earthquakes are used to quantify both compression (Q and shear(Q ) wave attenuation within the Lau backarc basin. A spectral-ratio method is employed to measure differential attenuation (t*) between the depth phases sS, pP, and sP and the direct S and P phases over the frequency band 0.05 and 0.5 Hz. We use a stacking algorithm to combine the spectra of several phase pairs from a single event, having similar azimuth and range, to obtain more robust t* measurements; these estimates are then used to compute the average Q above the focal depth. Q and Q are measured directly from the sS-S and pP-P phase pairs respectively, however, the interpretation of t* measured from sP-P requires assumptions about the ratio Q ?/Q . We find an empirical ratio of Q ?/Q = 1.93 for this region and use it to compute Q and Q from the Q _sP observations. We observe lateral and depth variations in both Q _β and Q _>alpha< beneath the tectonically active Lau Basin and the geologically older, inactive Lau Ridge and Fiji Plateau. The upper 200 km beneath the Central and Northern Lau Basin show a Q of 45–57 and a Q of 102–121, and Q appears to increase rapidly with depth. The upper 600 km beneath the Lau backarc basin has a Q of 118–138, while over the same depth interval we observe a higher Q of 139–161 beneath the Lau Ridge and Fiji Plateau. We also find Q of 235–303 beneath the northern Lau Basin and a higher Q of 292–316 beneath the Fiji Plateau and the Lau Ridge measured directly from pP-P phase pairs. These geographic trends in the broadband Q measurements correlate with our previous long-period estimates of Q in this region, however, the broadband measurements themselves are higher by about a factor of two. These observations suggest substantial frequency dependence of Q in the upper mantle, beginning at frequencies less than 1.0 Hz and consistent with the power-law form Q∝ω^α with between -0.1 and -0.3.     

<Emphasis Type="Italic">Q</Emphasis><Subscript>c</Subscript> and <Emphasis Type="Italic">Q</Emphasis><Subscript>S</Subscript> wave attenuation of South African earthquakes

Quality factor Q, which describes the attenuation of seismic waves with distance, was determined for South Africa using data recorded by the South African National Seismograph Network. Because of an objective paucity of seismicity in South Africa and modernisation of the seismograph network only in 2007, I carried out a coda wave decay analysis on only 13 tectonic earthquakes and 7 mine-related events for the magnitude range 3.6?≤?M _L ?≤?4.4. Up to five seismograph stations were utilised to determine Q _c for frequencies at 2, 4, 8 and 16 Hz resulting in 84 individual measurements. The constants Q ₀ and α were determined for the attenuation relation Q _c(f)?=?Q ₀ f ^α . The result was Q ₀?=?396?±?29 and α?=?0.72?±?0.04 for a lapse time of 1.9*(t _s???t ₀) (time from origin time t ₀ to the start of coda analysis window is 1.9 times the S-travel time, t _s) and a coda window length of 80 s. This lapse time and coda window length were found to fit the most individual frequencies for a signal-to-noise ratio of at least 3 and a minimum absolute correlation coefficient for the envelope of 0.5. For a positive correlation coefficient, the envelope amplitude increases with time and Q _c was not calculated. The derived Q _c was verified using the spectral ratio method on a smaller data set consisting of nine earthquakes and one mine-related event recorded by up to four seismograph stations. Since the spectral ratio method requires absolute amplitudes in its calculations, site response tests were performed to select four appropriate stations without soil amplification and/or signal distortion. The result obtained for Q _S was Q ₀?=?391?±?130 and α?=?0.60?±?0.16, which agrees well with the coda Q _c result.     

A numerical study on the occurrence of flowing wells in the discharge area of basins due to the upward hydraulic gradient induced wellbore flow

The occurrence of flowing wells in basins has been found to be closely related to the discharge area with an upward hydraulic gradient. Unfortunately, previous studies on upward gradient induced wellbore flow with equaling total inflow (Q_in) in the deep and total outflow (Q_out) in the shallow could not explain the occurrence of flowing wells. By representing wells using the MNW2 Package imbedded in MODFLOW 2005, we obtain the exchange of groundwater between the aquifer and the well in the discharge area of 3D unit basins and identify three scenarios: Q_in = Q_out, Q_in > Q_out > 0 and Q_in > Q_out = 0. The relationship of Q_in > Q_out well explains why flowing wells only develop in a limited part of the discharge area. Sensitivity analysis shows that well location, water table undulation, and basin length–depth ratio do not change the profile of the ratio of cumulative flow rate in a flowing well to total inflow (Q_v/Q_in) versus the relative elevation in the inflow segment, z_in*, but could significantly change the length of the inflow segment; well depth could change both the length of the inflow segment and the profile of Q_v/Q_in versus z_in*. Based on the numerical results in homogeneous and isotropic basins with different basin length–depth ratios in the current study, the ratio of inflow in the lower half part of a flowing well to the total inflow is found to be at least 67% and could be close to 100%, indicating that water at the outlets of flowing wells with long open sections is mainly from the deep part of the well.