SWAT hydrological model as a DaaS cloud service

Earth Science community depends on the exploration, analysis and reprocessing of high volumes of data as well as the modeling and simulation of complex coupled systems on multiple scales. The main aim of this article is to introduce a new hydrological modeling service based on the Soil and Water Assessment Tool (SWAT) (Arnold et al. J American Water Resour Assoc 34(1), 73–89, 1998 ; Arnold and Fohrer Hydrol Process 19(3), 563–572, 2005) model using high efficiency, resource sharing and low cost cloud computing resources (Astsatryan et al. International Journal of Scientific & Engineering Research 1(1), 1130–1133, 2014). Such a Desktop as a Service (DaaS) approach allowing users to work from anywhere, and gives centralized desktop management and great performance. Within the Spatial Data Infrastructure (SDI) and cloud platform, the DaaS service gives secure access to the model and a centralized data storage to get a SWAT model input. The article illustrates the analyses of the implementation of the SWAT model for the Sotk watershed of Lake Sevan in Armenia (Sargsyan 2007).     

An enhanced SMA based SCS-CN inspired model for watershed runoff prediction

Incorporation of initial soil moisture (V ₀) in the Soil Conservation Service Curve Number (SCS-CN) methodology helps to avoid the sudden jumps in Curve Number (CN) and, in turn, in computed runoff. It invoked the development of an enhanced (yet simple) Soil Moisture Accounting (SMA) procedure-based-SCS-CN inspired model, by incorporating initial moisture (V ₀). Its performance is tested using a dataset of 152 small to large watersheds of USDA (total 38,169 storm events), and compared with original SCS-CN method, Mishra and Singh (Acta Geophys Polon 50(3):457–477, 2002), Michel et al. (Water Resour Res 41(2):W02011, 2005) and Singh et al. (Water Resour Manag 29(11): 4111–4127, 2015) model using four statistical indices (RMSE, R ², PBIAS and NSE) and rank grading system (RGS). The proposed model scores highest (= 691 marks out of maximum 2280 marks) (Rank I) followed by Singh et al. (Water Resour Manag 29(11):4111–4127, 2015) model with 642 marks (Rank II), Michel et al. (Water Resour Res 41(2):W02011, 2005) model with 376 marks (Rank III) and Mishra and Singh model with 362 marks (= Rank IV). The original SCS-CN model, however, performs the poorest of all with 209 marks (Rank V).     

The vulnerability of the European air traffic network to spatial hazards

The 2010 eruption of the Eyjafjallajökull volcano had a devastating effect on the European air traffic network, preventing air travel throughout most of Europe for 6 days (Oroian in ProEnvironment 3:5–8, 2010). The severity of the disruption was surprising as previous research suggests that this type of network should be tolerant to random hazard (Albert et al. in Nature 406(6794):378–382, 2000; Strogatz in Nature 410(6825):268–276, 2001). The source of this hazard tolerance lies in the degree distribution of the network which, for many real-world networks, has been shown to follow a power law (Albert et al. in Nature 401(6749):130–131, 1999; Albert et al. in Nature 406(6794):378–382, 2000). In this paper, we demonstrate that the ash cloud was unexpectedly disruptive because it was spatially coherent rather than uniformly random. We analyse the spatial dependence in air traffic networks and demonstrate how the combination of their geographical distribution and their network architectures jeopardises their inherent hazard tolerance.     

Qanat: a resource or a hazard? Rebuttal to “Qanat is not a hazard” by Parise (Environ Earth Sci 2016 75: 1476)

A reply essay is presented on the rebuttal article by Parise (Environ Earth Sci 75(23):1476, 2016) suggesting that qanat is not a hazard. It is presented as a refutation on the paper by Abbasnejad et al. (Environ Earth Sci 75:1306, 2016) in which the authors have explained the etiology, the impacts and remedies of a qanat hazard in Iranian urban areas. Since the majority of qanats in Iranian urban areas are abandoned, useless and threatening, according to definition and in comparison with similar features introduced as a hazard, they are considered as a hazard too. However, this does not mean, and Abbasnejad et al. (in Environ Earth Sci 75:1306, 2016) have not claimed, that all qanats are hazardous. In addition, the authors who have studied qanats in Iranian urban areas, before Abbasnejad et al. (in Environ Earth Sci 75:1306, 2016), have also considered qanats as a hazard.     

Discussion of “From saturated to unsaturated conditions and vice versa” by Martí Lloret-Cabot,Simon J. Wheeler,Jubert A. Pineda,Enrique Romero,and Daichao Sheng (https://doi.org/10.1007/s11440-017-0577-6)

The authors (Lloret-Cabot et al. in Acta Geotech 1–23, 2017) applied the glasgow coupled model (GCM), originally proposed by Wheeler et al. (Géotechnique 53(1):41–54, 2003), to the simulation of several experimental tests that involve transition between saturated and unsaturated states. The authors show qualitatively, but not quantitatively, predictions of GCM for shrinkage during air drying of normally consolidated samples (Fig. 13) under low mechanical stress conditions, without presenting the material parameters. The discussers, who have worked with GCM to model multilayer deposition of tailings/soft soils (Qi in Numerical investigation for slope stability of expansive soils and large strain consolidation of soft soils. Doctoral dissertation, University of Ottawa, 2017; Qi et al. in J Geotech Geoenviron Eng 143(7):04017018, 2017, Qi et al. J Geotech Geoenviron Eng 143(7):04017019, 2017), have made quantitative predictions of similar cases. Satisfactory simulations of such cases using GCM are sensitive to the selection of the coupling parameters k₁ and k₂. By considering an alternative analytical form of GCM, an analytical procedure can be derived for calibrating the coupling parameters for problems involving virgin drying.     

The formation of porphyry copper deposits

Copper is a moderately incompatible chalcophile element. Its behavior is strongly controlled by sulfides. The speciation of sulfur is controlled by oxygen fugacity. Therefore, porphyry Cu deposits are usually oxidized (with oxygen fugacities > ΔFMQ +2) (Mungall 2002; Sun et al. 2015). The problem is that while most of the magmas at convergent margins are highly oxidized, porphyry Cu deposits are very rare, suggesting that high oxygen fugacity alone is not sufficient. Partial melting of mantle peridotite even at very high oxygen fugacities forms arc magmas with initial Cu contents too low to form porphyry Cu deposits directly (Lee et al. 2012; Wilkinson 2013). Here we show that partial melting of subducted young oceanic slabs at high oxygen fugacity (>ΔFMQ +2) may form magmas with initial Cu contents up to >500 ppm, favorable for porphyry mineralization. Pre-enrichment of Cu through sulfide saturation and accumulation is not necessarily beneficial to porphyry Cu mineralization. In contrast, re-melting of porphyritic hydrothermal sulfide associated with iron oxides may have major contributions to porphyry deposits. Thick overriding continental crust reduces the “leakage” of hydrothermal fluids, thereby promoting porphyry mineralization. Nevertheless, it is also more difficult for ore forming fluids to penetrate the thick continental crust to reach the depths of 2–4 km where porphyry deposits form.     

Micro-mechanical failure analysis of wet granular matter

We employ a novel fluid–particle model to study the shearing behavior of granular soils under different saturation levels, ranging from the dry material via the capillary bridge regime to higher saturation levels with percolating clusters. The full complexity of possible liquid morphologies Scheel et al. (Nat Mater 7(3):189–193, 2008. doi: 10.1038/nmat2117) is taken into account, implying the formation of isolated arbitrary-sized liquid clusters with individual Laplace pressures that evolve by liquid exchange via films on the grain surface Melnikov et al. (Phys Rev E 92(022):206, 2015. doi: 10.1103/PhysRevE.92.022206). Liquid clusters can grow in size, shrink, merge and split, depending on local conditions, changes of accessible liquid and the pore space morphology determined by the granular phase. This phase is represented by a discrete particle model based on contact dynamics Brendel et al. (Contact dynamics for beginners. Wiley-VCH, Weinheim, 2005. doi: 10.1002/352760362X.ch14), where capillary forces exerted from a liquid phase add to the motion of spherical particles. We study the macroscopic response of the system due to an external compression force at various liquid contents with the help of triaxial shear tests. Additionally, the change in liquid cluster distributions during the compression due to the deformation of the pore space is evaluated close to the critical load.     

Garnet–orthopyroxene (GOPX) geothermometer: a comparative study

The garnet–orthopyroxene pairs are commonly found in the assemblages of basic granulites/charnockite and hence are suitable for estimating equilibrium temperature of these metamorphic rocks. At present, there are many calibrations of garnet–orthopyroxene thermometer that may confuse geologists in choosing a reliable thermometer. To test the accuracy of the garnet–orthopyroxene thermometers, we have applied 14 models formulated by a number of workers since 1980 to date. We have collated 51 samples from the literature all over the world, which has been processed through the “Gt-Opx.EXE” software. Based on the present study, we have identified a set of the best among all the 14 models which were considered under this comparative study. We have concluded that the five garnet–orthopyroxene (Gt-Opx) thermometers are the most valid and reliable of this kind of thermometer: Aranovich and Berman (Am Mineral 82:345–353, 1997), Raith et al. (Earth Sci 73:211–244, 1983), Harley (Contrib Mineral Petrol 86:359–373, 1984), Nimis and Grütter (Contrib Mineral Petrol 159:411–427, 2010), and Sen and Bhattacharya (Contrib Mineral Petrol 88:64–71, 1984).     

Comparative approach of three popular intrinsic vulnerability methods: case of the Beni Amir groundwater (Morocco)

The study of the intrinsic vulnerability of groundwater resources to pollution is an effective tool to control their quality degradation and contribute to their protection. It is used to delimit the vulnerable zones which do not withstand a large flow of pollutants introduced from the soil surface. Three methods of assessing the intrinsic vulnerability of groundwater: DRASTIC (Depth to water table, Recharge, Aquifer, Soil type, Topography, Impact of zone vadose, Hydraulic conductivity), DRSTI, and GOD (Groundwater occurrence, Overall aquifer class and Depth of water table) coupled with a geographic information system (GIS) are applied to the groundwater of Beni Amir, and they are compared in order to adopt the method which better characterizes the vulnerability of the aquifer to pollution. The validation of this application was made by measurements of the nitrate levels in the aquifer. Because the pollution of groundwater, in this plain, is a direct consequence of agricultural activities characterized by an intensive fertilizer application. The results clearly show that the rate of the coincidence, between the measured nitrate concentrations and the different classes of vulnerability of three methods, is 81.81, 54.54, 72.72, and 27.27%, respectively, for methods DRASTIC (classification of Engel et al. (Int Assoc Hydrol Sci Publi 235:521–526, 1996)), DRASTIC (classifications of Aller et al. (1987)), DRSTI, and GOD. Of this rate of coincidences, the DRASTIC method, with the classification of Engel et al. (Int Assoc Hydrol Sci Publi 235:521–526, 1996), allows a finer assessment and turns out the most representative of the study area.     

Discussion: “Xenoliths in ultrapotassic volcanic rocks in the Lhasa block: direct evidence for crust–mantle mixing and metamorphism in the deep crust” by Wang et al. 2016 (Contributions to Mineralogy and Petrology) 171:62

Wang et al. (Contrib Mineral Petrol 171:62, 2016a) present data on composition of xenolith from Southern Tibet and conclude that ulrapotassic melts from the region formed by melting mantle, and complex interaction with a crustal component. In this discussion we demonstrate that numerous observations presented by Wang et al. (2016a) can be explained by partial melting of crust followed by interaction between that melt and the mantle. We show that this model can explain the variability of magmas in such suits without evoking occurrence of coincidental, unrelated events. Moreover we demonstrate that our model of a crustal origin of the proto-shoshonite melts is now supported by independent lines of evidence such as geochemistry of restites after high- and ultrahigh- pressure melting and melt inclusion studies.     

Anisotropic mechanical properties of zircon and the effect of radiation damage

This study provides new insights into the relationship between radiation-dose-dependent structural damage due to natural U and Th impurities and the anisotropic mechanical properties (Poisson’s ratio, elastic modulus and hardness) of zircon. Natural zircon samples from Sri Lanka (see Muarakami et al. in Am Mineral 76:1510–1532, 1991) and synthetic samples, covering a dose range of zero up to 6.8 × 10¹⁸ α-decays/g, have been studied by nanoindentation. Measurements along the [100] crystallographic direction and calculations, based on elastic stiffness constants determined by Özkan (J Appl Phys 47:4772–4779, 1976), revealed a general radiation-induced decrease in stiffness (~54 %) and hardness (~48 %) and an increase in the Poisson’s ratio (~54 %) with increasing dose. Additional indentations on selected samples along the [001] allowed one to follow the amorphization process to the point that the mechanical properties are isotropic. This work shows that the radiation-dose-dependent changes of the mechanical properties of zircon can be directly correlated with the amorphous fraction as determined by previous investigations with local and global probes (Ríos et al. in J Phys Condens Matter 12:2401–2412, 2000a; Farnan and Salje in J Appl Phys 89:2084–2090, 2001; Zhang and Salje in J Phys Condens Matter 13:3057–3071, 2001). The excellent agreement, revealed by the different methods, indicates a large influence of structural and even local phenomena on the macroscopic mechanical properties. Therefore, this study indicates the importance of acquiring better knowledge about the mechanical long-term stability of radiation-damaged materials.     

Simplified SMA-inspired 1-parameter SCS-CN model for runoff estimation

This study proposes a simplified 1-parameter SCS-CN model (M5) based on Mishra-Singh (2002) model and soil moisture accounting (SMA) procedure for surface runoff estimation and compares its performance with the existing SCS-CN method (SCS, 1956) (M1), Michel 1-P model (Water Resour Res 41:1-6, 2005) (M2), Sahu 1-P model (Hydrol Process 21:2872-2881, 2007) (M3), and Ajmal et al. model (J Hydrol 530:623-633, 2015) (M4) using large rainfall–runoff dataset of 48,763 events from123 USDA-ARS watersheds. The performance of models was evaluated using three statistical error indices such as Nash-Sutcliffe efficiency (NSE), root mean square error (RMSE), percentage bias (PBIAS), and rank and grading system (RGS). Based on the results obtained, the models can be ranked as follows: M5?>?M4?>?M3?>?M1?>?M2, i.e., model M5 outperformed all the remaining four models M1–M4 and hence is recommended for field applications.     

Machine learning for graph-based representations of three-dimensional discrete fracture networks

Structural and topological information play a key role in modeling flow and transport through fractured rock in the subsurface. Discrete fracture network (DFN) computational suites such as dfnWorks (Hyman et al. Comput. Geosci. 84, 10–19 2015) are designed to simulate flow and transport in such porous media. Flow and transport calculations reveal that a small backbone of fractures exists, where most flow and transport occurs. Restricting the flowing fracture network to this backbone provides a significant reduction in the network’s effective size. However, the particle-tracking simulations needed to determine this reduction are computationally intensive. Such methods may be impractical for large systems or for robust uncertainty quantification of fracture networks, where thousands of forward simulations are needed to bound system behavior. In this paper, we develop an alternative network reduction approach to characterizing transport in DFNs, by combining graph theoretical and machine learning methods. We consider a graph representation where nodes signify fractures and edges denote their intersections. Using random forest and support vector machines, we rapidly identify a subnetwork that captures the flow patterns of the full DFN, based primarily on node centrality features in the graph. Our supervised learning techniques train on particle-tracking backbone paths found by dfnWorks, but run in negligible time compared to those simulations. We find that our predictions can reduce the network to approximately 20% of its original size, while still generating breakthrough curves consistent with those of the original network.     

Numerical analysis of effect of baffle configuration on impact force exerted from rock avalanches

In mountainous areas, channelized rock avalanches swarm downslope leading to large impact forces on building structures in residential areas. Arrays of rock avalanche baffles are usually installed in front of rigid barriers to attenuate the flow energy of rock avalanches. However, previous studies have not sufficiently addressed the mechanisms of interaction between the rock avalanches and baffles. In addition, empirical design approaches such as debris flow (Tang et al., Quat Int 250:63–73, 2012), rockfall (Spang and Rautenstrauch, 1237–1243, 1988), snow avalanches (Favier et al., 14:3–15, 2012), and rock avalanches (Manzella and Labiouse, Landslides 10:23–36, 2013), which are applied in natural geo-disasters mitigation cannot met construction requirements. This study presents details of numerical modeling using the discrete element method (DEM) to investigate the effect of the configuration of baffles (number and spacing of baffle columns and rows) on the impact force that rock avalanches exert on baffles. The numerical modeling is firstly conducted to provide insights into the flow interaction between rock avalanches and an array of baffles. Then, a modeling analysis is made to investigate the change pattern of the impact force with respect to baffle configurations. The results demonstrate that three crucial influencing factors (baffle row numbers, baffle column spacing, and baffle row spacing) have close relationship with energy dissipation of baffles. Interestingly, it is found that capacity of energy dissipation of baffles increases with increasing baffle row numbers and baffle row spacing, while it decreases with increasing baffle column spacing. The results obtained from this study are useful for facilitating design of baffles against rock avalanches.     

Reply to the comments on “Xenoliths in ultrapotassic volcanic rocks in the Lhasa block: direct evidence for crust–mantle mixing and metamorphism in the deep crust”

Stepanov et al. (Contrib Mineral Petrol, 2017) question our conclusion that the UPVs in southern Tibet were derived by partial melting of an old, metasomatized subcontinental lithospheric mantle (SCLM) of the subducted Indian plate. Instead, they propose that these ultrapotassic volcanic rocks (UPVs) are shoshonitic and were generated in two steps: direct melting of crustal rocks first, and then the melts interacted with mantle peridotite. However, the trace element, isotopic, thermal, structural, and seismic evidence is consistent with the xenolith evidence (Wang et al in Contrib Mineral Petrol 172:62, 2016) for hybridisation of ascending Indian subcontinental lithospheric mantle-derived UPV magmas with the deep, isotopically unevolved, Tibetan crust. This necessitates a model whereby partial melting of subducting Indian SCLM generates the UPV suite of southern Tibet.     

Shock trains on a planar beach: quasi-analytical and fully numerical solutions

This study, part of the Special Issue dedicated to the 70th anniversary of Professor Efim Pelinovsky, focuses on a topic that has been central in Professor Pelinovsky’s research, i.e. the analytical and numerical modelling of shallow water waves. We specifically focus on the evolution of trains of shock waves on a planar beach. Antuono (J Fluid Mech 658:166–187, 2011) has, for the first time, proposed a quasi-analytical solution for a train of shock waves forced by a constant Riemann invariant. The present contribution clarifies the validity of such solution and its value for benchmarking nonlinear shallow water equation solvers. Hence, the same tests of Antuono (J Fluid Mech 658:166–187, 2011) have been run by means of the solver of Brocchini et al. (Coast Eng 43(2):105–129, 2001) revealing surprisingly and reassuring good agreements. This provides significant support to the mentioned analytical solution and allows to critically analyse the eventual discrepancies, due to the practicalities of running numerical shallow water solutions (e.g. influence of the boundary conditions, of the numerical resolution, etc.).     

Indirect unstructured hex-dominant mesh generation using tetrahedra recombination

Corner-point gridding is widely used in reservoir and basin modeling but generally yields approximations in the representation of geological interfaces. This paper introduces an indirect method to generate a hex-dominant mesh conformal to 3D geological surfaces and well paths suitable for finite-element and control-volume finite-element simulations. By indirect, we mean that the method first generates an unstructured tetrahedral mesh whose tetrahedra are then merged into primitives (hexahedra, prisms, and pyramids). More specifically, we focus on determining the optimal set of primitives that can be recombined from a given tetrahedral mesh. First, we detect in the tetrahedral mesh all the feasible volumetric primitives using a pattern-matching algorithm (Meshkat and Talmor Int. J. Numer. Meth. Eng. 49(1-2), 17–30 2000) that we re-visit and extend with configurations that account for degenerated tetrahedra (slivers). Then, we observe that selecting the optimal set of primitives among the feasible ones can be formalized as a maximum weighted independent set problem (Bomze et al. 1999), known to be \(\mathcal {N}\mathcal {P}\)-Complete. We propose several heuristic optimizations to find a reasonable set of primitives in a practical time. All the tetrahedra of each selected primitive are then merged to build the final unstructured hex-dominant mesh. This method is demonstrated on 3D geological models including a faulted and folded model and a discrete fracture network.     

Miocene echinoids from the Sadat area,south Gebel Ataqa,NW Gulf of Suez,Egypt: systematics,palaeobiogeography

The echinoid fauna from the Miocene sedimentary succession cropping out south Wadi Tweirig, and Wadi Hommath, south Gebel Ataqa, NW Gulf of Suez, has been examined with the aim to known their stratigraphic and paleogeographic distribution. The Miocene succession includes two formations: Sadat Formation, unconformably overlying the middle/upper Eocene rocks at the base and Hommath Formation at the top. Twenty-eight echinoid species (8 regular and 20 irregular) belonging to 18 genera, 13 families, and 7 orders have been identified, systematically described, and illustrated in this work. Eleven species are recorded for the first time from Egypt: ten of these came from the Hommath Formation (Schizechinus cf. serresii Desor (1856), Schizechinus pentagonus Kier 1972, Clypeaster cf. martini des Moulins 1837, Scutella checchiae occidentalis Desio 1934, Scutella melitensis Airaghi 1902, Echinodiscus desori Duncan and Sladen 1883, Echinolampas cf. zeitensis Fourtau 1920, Schizaster lovisatoi Cotteau 1895, Agassizia (Agassizia) powersi Kier 1972, and Hemipatagus ocellatus Defrance (1827)), and one from the Sadat Formation (Clypeaster campanulatus Schlotheim (1820)). The identified fauna shows a strong affinity with the Mediterranean bio-province.     

Ultimate Bearing Capacity Prediction of Eccentrically Inclined Loaded Strip Footings

Extensive laboratory model tests have been carried out on a strip footing resting over dry sand bed subjected to eccentrically inclined load to determine the ultimate bearing capacity (Patra et al. in Int J Geotech Eng 6(3):343–352, 2012a.  https://doi.org/10.3328/IJGE.2012.06.03.343-352, Int J Geotech Eng 6(4):507–514, b.  https://doi.org/10.3328/IJGE.2012.06.04.507-514). Similarly, lower bound calculations based on finite element method were performed to compute the bearing capacity of a strip footing subjected to an eccentric and inclined load lying over a cohesionless soil with varying embedment depth and relative density (Krabbenhoft et al. in Int J Geomech ASCE, 2014.  https://doi.org/10.1061/(ASCE)GM.1943-5622.0000332). The load may be applied in two ways namely, towards the center line and away from the center line of the footing. Based on the results (both experimental and numerical analyses), a neural network model is developed to predict the reduction factor that will be used in computing the ultimate bearing capacity of an eccentrically inclined loaded strip footing. This reduction factor (RF) is the ratio of the ultimate bearing capacity of the footing subjected to an eccentrically inclined load to the ultimate bearing capacity of the footing subjected to a centric vertical load. A thorough sensitivity analysis is carried out to evaluate the parameters affecting the reduction factor. Based on the weights of the developed neural network model, a neural interpretation diagram is developed to find out whether the input parameters have direct or inverse effect on the output. A prediction model equation is framed with the trained weights of the neural network as the model parameters. The predictions from ANN, and those from other approaches, are compared with the results computed from both experimentation and FEM analyses. The ANN model results are found to be more accurate and well matched with other results.     

Hybrid discretization of multi-phase flow in porous media in the presence of viscous,gravitational, and capillary forces

Multi-phase flow in porous media in the presence of viscous, gravitational, and capillary forces is described by advection diffusion equations with nonlinear parameters of relative permeability and capillary pressures. The conventional numerical method employs a fully implicit finite volume formulation. The phase-potential-based upwind direction is commonly used in computing the transport terms between two adjacent cells. The numerical method, however, often experiences non-convergence in a nonlinear iterative solution due to the discontinuity of transmissibilities, especially in transition between co-current and counter-current flows. Recently, Lee et al. (Adv. Wat. Res. 82, 27–38, 2015) proposed a hybrid upwinding method for the two-phase transport equation that comprises viscous and gravitational fluxes. The viscous part is a co-current flow with a one-point upwinding based on the total velocity and the buoyancy part is modeled by a counter-current flow with zero total velocity. The hybrid scheme yields C¹-continuous discretization for the transport equation and improves numerical convergence in the Newton nonlinear solver. Lee and Efendiev (Adv. Wat. Res. 96, 209–224, 2016) extended the hybrid upwind method to three-phase flow in the presence of gravity. In this paper, we present the hybrid-upwind formula in a generalized form that describes two- and three-phase flows with viscous, gravity, and capillary forces. In the derivation of the hybrid scheme for capillarity, we note that there is a strong similarity in mathematical formulation between gravity and capillarity. We thus greatly utilize the previous derivation of the hybrid upwind scheme for gravitational force in deriving that for capillary force. Furthermore, we also discuss some mathematical issues related to heterogeneous capillary domains and propose a simple discretization model by adapting multi-valued capillary pressures at the end points of capillary pressure curves. We demonstrate this new model always admits a consistent solution that is within the discretization error. This new generalized hybrid scheme yields a discretization method that improves numerical stability in reservoir simulation.