Large Viscosity Variations Research Articles

The Influence of Large Variations in Fluid Density and Viscosity on the Resonance Characteristics of Tuning Forks Simulated by Finite Element Method

The use of tuning forks to measure fluid density and viscosity is widely employed in fields such as food, medicine, textiles, automobiles, petrochemicals, and deep drilling. The explicit analytical model based on the Euler–Bernoulli cantilever-beam theory for the relationship between tuning-fork resonance characteristics and the density and viscosity of fluid is only applicable to the situation where the fluid viscous effect is very small. In this paper, the finite element method is used to simulate the influence of large variations in fluid density and viscosity on the resonance characteristic parameters (resonant frequency and quality factor) of the tuning fork. The numerical simulation results are compared with the analytical analysis results and experimental measurement results. Then, the sensitivity of tuning-fork resonance characteristic parameters to fluid density and viscosity is studied. The results show that compared with the analytical results, the numerical simulation results have a higher degree of agreement with the experimental measurement results. The relative difference in resonant frequency is less than 2%, and the relative difference in quality factor is less than 4%. This indicates that the finite element method includes the influence of fluid viscosity on tuning-fork resonance parameters, which is more in line with the actual conditions than the analytical model. Simulating and analyzing the sensitivity of the tuning fork to fluid density and viscosity by the finite element method, it is possible to consider the situation where fluid density and viscosity vary over a large range. Compared with experimental measurements, this method has higher efficiency and can significantly save time and economic costs. This study can overcome the limitation of existing explicit analytical models, which are only applicable when the viscous effects of the fluid are very small. It enables a more accurate simulation of the coupling vibration between tuning forks and fluids, thereby providing theoretical references for further optimizing tuning-fork structural parameters to enhance the accuracy of measuring fluid characteristic parameters.

Modeling compressed turbulent plasma with rapid viscosity variations

We propose two-equation models in order to capture the dynamics of a turbulent plasma undergoing compression and experiencing large viscosity variations. The models account for possible relaminarization phases and rapid viscosity changes through closures dependent on the turbulent Reynolds and on the viscosity Froude numbers. These closures are determined from a data-driven approach using eddy-damped quasi-normal Markovian simulations. The best model is able to mimic the various self-similar regimes identified in Viciconte et al. [Phys. Rev. E 97, 023201 (1998)] and to recover the rapid transition limits identified by G. N. Coleman and N. N. Mansour [Phys. Fluids A 3, 2255 (1991)].

Evolutionary divergence of anaphase spindle mechanics in nematode embryos constrained by antagonistic pulling and viscous forces.

Cellular functions such as cell division are remarkably conserved across phyla. However, the evolutionary principles of cellular organization that drive them are less well explored. Thus, an essential question remains: to what extent do cellular parameters evolve without altering the basic functions they sustain? Here we have observed six different nematode species for which the mitotic spindle is positioned asymmetrically during the first embryonic division. Whereas the C. elegans spindle undergoes oscillations during its displacement, the spindle elongates without oscillations in other species. We asked which evolutionary changes in biophysical parameters could explain differences in spindle motion while maintaining a constant output. Using laser microsurgery of the spindle, we revealed that all species are subjected to cortical pulling forces of varying magnitudes. Using a viscoelastic model to fit the recoil trajectories and with an independent measurement of cytoplasmic viscosity, we extracted the values of cytoplasmic drag, cortical pulling forces, and spindle elasticity for all species. We found large variations in cytoplasmic viscosity, whereas cortical pulling forces and elasticity were often more constrained. In agreement with previous simulations, we found that increased viscosity correlates with decreased oscillation speeds across species. However, the absence of oscillations in some species despite low viscosity can only be explained by smaller pulling forces. Consequently, we find that spindle mobility across the species analyzed here is characterized by a tradeoff between cytoplasmic viscosity and pulling forces normalized by the size of the embryo. Our work provides a framework for understanding mechanical constraints on evolutionary diversification of spindle mobility.

Mathematical modelling of mantle convection at a high Rayleigh number with variable viscosity and viscous dissipation

In this paper, the classical Rayleigh–Bénard convection model is considered and solved numerically for extremely large viscosity variations (i.e., up to 10^{30}) across the mantle at a high Rayleigh number. The Arrhenius form of viscosity is defined as a cut-off viscosity function. The effects of viscosity variation and viscous dissipation on convection with temperature-dependent viscosity and also temperature- and pressure-dependent viscosity are shown through the figures of temperature profiles and streamline contours. The values of Nusselt number and root mean square velocity indicate that the convection becomes significantly weak as viscosity variation and viscous dissipation are increased at a fixed pressure dependence parameter.

The effect of lateral variations in Earth structure on Last Interglacial sea level

SUMMARYIt is generally agreed that the Last Interglacial (LIG; ∼130–115 ka) was a time when global average temperatures and global mean sea level were higher than they are today. However, the exact timing, magnitude and spatial pattern of ice melt is much debated. One difficulty in extracting past global mean sea level from local observations is that their elevations need to be corrected for glacial isostatic adjustment (GIA), which requires knowledge of Earth’s internal viscoelastic structure. While this structure is generally assumed to be radially symmetric, evidence from seismology, geodynamics and mineral physics indicates that large lateral variations in viscosity exist within the mantle. In this study, we construct a new model of Earth’s internal structure by converting shear wave speed into viscosity using parametrizations from mineral physics experiments and geodynamic constraints on Earth’s thermal structure. We use this 3-D Earth structure, which includes both variations in lithospheric thickness and lateral variations in viscosity, to calculate the first 3-D GIA prediction for LIG sea level. We find that the difference between predictions with and without lateral Earth structure can be metres to 10s of metres in the near field of former ice sheets, and up to a few metres in their far field. We demonstrate how forebulge dynamics and continental levering are affected by laterally varying Earth structure, with a particular focus on those sites with prominent LIG sea level records. Results from four 3-D GIA calculations show that accounting for lateral structure can act to increase local sea level by up to ∼1.5 m at the Seychelles and minimally decrease it in Western Australia. We acknowledge that this result is only based on a few simulations, but if robust, this shift brings estimates of global mean sea level from these two sites into closer agreement with each other. We further demonstrate that simulations with a suitable radial viscosity profile can be used to locally approximate the 3-D GIA result, but that these radial profiles cannot be found by simply averaging viscosity below the sea level indicator site.

Linear superposition time for variable-rate, variable-pressure production from a fractured tight gas well

Hydraulically fractured, tight gas wells produce under large pressure drawdowns, which may cause three-to ten-fold variations in the gas viscosity-compressibility product, particularly in the vicinity of the fracture, over the life of the production. Under such a large variation, the pseudo-pressure transformation cannot remove nonlinearity of the gas diffusion equation. The remaining nonlinearity makes the application of superposition principal questionable; specifically the linear superposition time based on linear flow, a dominant flow regime for hydraulically fractured tight gas reservoirs. Not accounting for the effect of large gas viscosity-compressibility product variations in pressure and rate transient analysis of fractured, tight gas well performances leads to misinterpretation of reservoir characteristics.In this paper, a perturbation-Green's function solution is developed for the nonlinear gas diffusion equation. Because each term in the perturbation solution represents the solution for a linearized problem, term-by-term application of the superposition principal is permitted. The semi-analytical nature of the solution also enables us to derive approximations for practical use. This paper explains the solution procedure, compares and verifies the solution with existing simulators presents variable-rate approximations in terms of new superposition time and discusses the results to provide guidelines for the analysis of fractured, tight gas well performances under large variations of viscosity and compressibility.

Largest WAG Pilot in Giant Al-Shaheen Field Reveals Optimization Methods

This article, written by JPT Technology Editor Chris Carpenter, contains highlights of paper SPE 190343, “Challenges and Learnings From Operating the Largest Offshore WAG in the Giant Al-Shaheen Field and Ways To Optimize Future WAG Developments,” by M. Pal, S. Furqan Gilani, and G. Tarsauliya, North Oil Company, prepared for the 2018 SPE EOR Conference at Oil and Gas West Asia, Muscat, Oman, 26–28 March. The paper has not been peer reviewed. This paper discusses the operation of the largest offshore high-pressure water-alternating-gas (WAG) injection pilot using hydrocarbon (HC) gas, with a gas-injection capacity of approximately 100 MMscf/D, in the giant Al-Shaheen field offshore Qatar. The paper also analyzes the current HC-WAG pilot and optimizes future development scenarios for extending HC-WAG injection to the field scale. Al-Shaheen WAG Background As part of the field-development plan, a comprehensive enhanced oil recovery (EOR) screening study was performed that explored the range of available EOR techniques. On the basis of these studies, WAG injection was identified as one of the more-promising EOR options for the field. Additional studies performed in 2006 and 2007 indicated a large full-field potential for WAG in Al-Shaheen, with an estimated recovery of hundreds of millions of barrels over and above water-flood recovery. Data analysis and intensive simulation work, including both a history match and several forecast sensitivities, were conducted to understand and to estimate WAG potential on Al-Shaheen reservoirs. The Al-Shaheen WAG process is largely an immiscible process with some multi contact miscibility effects between injectors and producers. Overall, WAG has been a success in the field in terms of seeing incremental oil gains on several patterns. Although physically successful, the project had been challenging operationally. In spite of incremental gains, a number of operational challenges were encountered during the execution of the WAG pilot. These challenges mostly were related to WAG-cycle conversions and operational efficiency of gas-compression systems. WAG Recovery Optimization WAG will be a major part of the development of the mature-flank and heavy-oil areas of the field. The areas proposed for WAG are undersaturated with gas at reservoir conditions. Incremental recovery from WAG could be affected by thermodynamic effects caused by large variations in permeability, viscosity, and saturation pressure in the field along with several other variables described in the complete paper. To analyze the effect of these optimization parameters on incremental oil recovery from WAG, a sector model with varying static and dynamic properties was built (Fig. 1). The properties of the sector model were based on the history-matched model of the neighboring areas. The well spacing was 700 ft, which lies in the middle of the well-spacing ranges of the ongoing WAG pilot wells (500 to 1,000 ft). The WAG sector model was initialized with a range of fixed oil API gravities. Different thermodynamic effects of gas injection have been examined on the sector model. The model is also simulated for different WAG cycle lengths (1/1, 3/3, 6/6, and 12/12), ratios (1:2, 2:1, and 1:3), and slug sizes (1–30 years of WAG). The effects of choking production wells and changes in injection gas composition were also tested on the sector model.

Theoretical and numerical analysis of the lattice kinetic scheme for complex-flow simulations.

The lattice kinetic scheme (LKS) is a modified version of the classical single relaxation time lattice Boltzmann method. Although used for many applications, especially when large variations in viscosity are involved, a thorough analysis of the scheme has not been provided yet. In the context of this work, the macroscopic behavior of this scheme is evaluated through the Chapman-Enskog analysis. It is shown that the additional degree of freedom provided in the scheme allows for an independent control of higher-order moments. These results are further corroborated by numerical simulations. The behavior of this numerical scheme is studied for selected external and internal flows to clarify the effect of the free parameter on the different moments of the distribution function. It is shown that it is more stable than SRT (single relaxation time) when confronted to fully periodic under-resolved simulations (especially for λ≈1). It can also help minimize the error coming from the viscosity-dependence of the wall position when combined with the bounce-back approach; although still present, viscosity-dependence of the wall position is reduced. Furthermore, as shown through the multiscale analysis, specific choices of the free parameter can cancel out the leading-order error. Overall, the LKS is shown to be a useful and efficient alternative to the SRT method for simulating numerically complex flows.

Asthenosphere Flow Modulated by Megathrust Earthquake Cycles

AbstractSubduction megathrusts develop the largest earthquakes, often close to large population centers. Understanding the dynamics of deformation at subduction zones is therefore important to better assess seismic hazards. Here I develop consistent earthquake cycle simulations that incorporate localized and distributed deformation based on laboratory‐derived constitutive laws by combining boundary and volume elements to represent the mechanical coupling between megathrust slip and solid‐state flow in the oceanic asthenosphere and in the mantle wedge. The model is simplified, in two dimensions, but may help the interpretation of geodetic data. Megathrust earthquakes and slow‐slip events modulate the strain rate in the upper mantle, leading to large variations of effective viscosity in space and time and a complex pattern of surface deformation. While fault slip and flow in the mantle wedge generate surface displacements in the same, that is, seaward, direction, the viscoelastic relaxation in the oceanic asthenosphere generates transient surface deformation in the opposite, that is, landward, direction above the rupture area of the mainshock. Aseismic deformation above the seismogenic zone may be challenging to record, but it may reveal important constraints about the rheology of the subducting plate.

Constraints on dynamic topography from asymmetric subsidence of the mid-ocean ridges

Stresses from mantle convection deflect Earth's surface vertically, producing dynamic topography that is important for continental dynamics and sea-level change but difficult to observe due to overprinting by isostatic topography. For long wavelengths (∼104 km), the amplitude of dynamic topography is particularly uncertain, with mantle flow models typically suggesting larger amplitudes (>1000 m) than direct observations. Here we develop a new constraint on the amplitude of long-wavelength dynamic topography by examining asymmetries in seafloor bathymetry across mid-ocean ridges. We compare bathymetric profiles across the Mid-Atlantic Ridge (MAR) and the East Pacific Rise (EPR) and we find that the South American flank of both ridges subsides faster than its opposing flank. This pattern is consistent with dynamic subsidence across South America, supported by downwelling in the lower mantle. To constrain the amplitude of dynamic topography, we compare bathymetric profiles across both ridges after correcting bathymetry for several different models of dynamic topography with varying amplitudes and spatial patterns. We find that long-wavelength dynamic topography with an amplitude of only ∼500 m explains the observed asymmetry of the MAR. A similar model can explain EPR asymmetry but is complicated by additional asymmetrical topography associated with tectonic, crustal thickness, and/or asthenospheric temperature asymmetries across the EPR. After removing 500 m of dynamic topography, both the MAR and EPR exhibit a slower seafloor subsidence rate (∼280–290 m/Myr1/2) than previously reported. Our finding of only ∼500 m of long-wavelength dynamic topography may indicate the importance of thermochemical convection and/or large viscosity variations for lower mantle dynamics.

A Hybrid Approach to Data Assimilation for Reconstructing the Evolution of Mantle Dynamics

AbstractQuantifying past mantle dynamic processes represents a major challenge in understanding the temporal evolution of the solid earth. Mantle convection modeling with data assimilation is one of the most powerful tools to investigate the dynamics of plate subduction and mantle convection. Although various data assimilation methods, both forward and inverse, have been created, these methods all have limitations in their capabilities to represent the real earth. Pure forward models tend to miss important mantle structures due to the incorrect initial condition and thus may lead to incorrect mantle evolution. In contrast, pure tomography‐based models cannot effectively resolve the fine slab structure and would fail to predict important subduction‐zone dynamic processes. Here we propose a hybrid data assimilation approach that combines the unique power of the sequential and adjoint algorithms, which can properly capture the detailed evolution of the downgoing slab and the tomographically constrained mantle structures, respectively. We apply this new method to reconstructing mantle dynamics below the western U.S. while considering large lateral viscosity variations. By comparing this result with those from several existing data assimilation methods, we demonstrate that the hybrid modeling approach recovers the realistic 4‐D mantle dynamics the best.

Size Distribution and Dispersion of Droplets Generated by Impingement of Breaking Waves on Oil Slicks

AbstractThis laboratory experimental study investigates the temporal evolution of the size distribution of subsurface oil droplets generated as breaking waves entrain oil slicks. The measurements are performed for varying wave energy, as well as large variations in oil viscosity and oil‐water interfacial tension, the latter achieved by premixing the oil with dispersant. In situ measurements using digital inline holography at two magnifications are applied for measuring the droplet sizes and Particle Image Velocimetry (PIV) for determining the temporal evolution of turbulence after wave breaking. All early (2–10 s) size distributions have two distinct size ranges with different slopes. For low dispersant to oil ratios (DOR), the transition between them could be predicted based on a turbulent Weber (We) number in the 2–4 range, suggesting that turbulence plays an important role. For smaller droplets, all the number size distributions have power of about −2.1, and for larger droplets, the power decreases well below −3. The measured steepening of the size distribution over time is predicted by a simple model involving buoyant rise and turbulence dispersion. Conversely, for DOR 1:100 and 1:25 oils, the diameter of slope transition decreases from ∼1 mm to 46 and 14 µm, respectively, much faster than the We‐based prediction, and the size distribution steepens with increasing DOR. Furthermore, the concentration of micron‐sized droplets of DOR 1:25 oil increases for the first 10 min after entrainment. These phenomena are presumably caused by the observed formation and breakup oil microthreads associated with tip streaming.

Eliminating volatile-induced surface porosity during resin transfer molding of a benzoxazine/epoxy blend

We studied the mechanism of volatile-induced surface porosity formation during the resin transfer molding (RTM) of aerospace composites using a blended benzoxazine/epoxy resin, and identified reduction strategies based on material and processing parameters. First, the influence of viscosity and pressure on resin volatilization were determined. Then, in situ data was collected during molding using a lab-scale RTM system for different cure cycles and catalyst concentrations. Finally, the surface quality of molded samples was evaluated. The results show that surface porosity occurs when cure shrinkage causes a sufficient decrease in cavity pressure prior to resin vitrification. The combination of thermal gradients and rapid gelation can generate large spatial variations in viscosity, rendering the coldest regions of a mold susceptible to porosity formation. However, material and cure cycle modifications can alter the resin cure kinetics, making it possible to delay the pressure drop until higher viscosities are attained to minimize porosity formation.

Numerical studies of thermal convection with temperature- and pressure-dependent viscosity at extreme viscosity contrasts

Motivated by convection of planetary mantles, we consider a mathematical model for Rayleigh-Bénard convection in a basally heated layer of a fluid whose viscosity depends strongly on temperature and pressure, defined in an Arrhenius form. The model is solved numerically for extremely large viscosity variations across a unit aspect ratio cell, and steady solutions for temperature, isotherms, and streamlines are obtained. To improve the efficiency of numerical computation, we introduce a modified viscosity law with a low temperature cutoff. We demonstrate that this simplification results in markedly improved numerical convergence without compromising accuracy. Continued numerical experiments suggest that narrow cells are preferred at extreme viscosity contrasts, and this conclusion is supported by a linear stability analysis.

Effects of aortic irregularities on blood flow

Anatomic aortic anomalies are seen in many medical conditions and are known to cause disturbances in blood flow. Turner syndrome (TS) is a genetic disorder occurring only in females where cardiovascular anomalies, particularly of the aorta, are frequently encountered. In this study, numerical simulations are applied to investigate the flow characteristics in four TS patient- related aortic arches (a normal geometry, dilatation, coarctation and elongation of the transverse aorta). The Quemada viscosity model was applied to account for the non-Newtonian behavior of blood. The blood is treated as a mixture consisting of water and red blood cells (RBC) where the RBCs are modeled as a convected scalar. The results show clear geometry effects where the flow structures and RBC distribution are significantly different between the aortas. Transitional flow is observed as a jet is formed due to a constriction in the descending aorta for the coarctation case. RBC dilution is found to vary between the aortas, influencing the WSS. Moreover, the local variations in RBC volume fraction may induce large viscosity variations, stressing the importance of accounting for the non-Newtonian effects.

Advantages of a conservative velocity interpolation (CVI) scheme for particle-in-cell methods with application in geodynamic modeling.

The particle‐in‐cell method is generally considered a flexible and robust method to model the geodynamic problems with chemical heterogeneity. However, velocity interpolation from grid points to particle locations is often performed without considering the divergence of the velocity field, which can lead to significant particle dispersion or clustering if those particles move through regions of strong velocity gradients. This may ultimately result in cells void of particles, which, if left untreated, may, in turn, lead to numerical inaccuracies. Here we apply a two‐dimensional conservative velocity interpolation (CVI) scheme to steady state and time‐dependent flow fields with strong velocity gradients (e.g., due to large local viscosity variation) and derive and apply the three‐dimensional equivalent. We show that the introduction of CVI significantly reduces the dispersion and clustering of particles in both steady state and time‐dependent flow problems and maintains a locally steady number of particles, without the need for ad hoc remedies such as very high initial particle densities or reseeding during the calculation. We illustrate that this method provides a significant improvement to particle distributions in common geodynamic modeling problems such as subduction zones or lithosphere‐asthenosphere boundary dynamics.

On the role of thermal heterogeneities on the rheology of MgO under conditions of the Earth’s lower mantle

The Earth’s mantle is characterised by large thermal heterogeneities associated with hot rising plumes and cold downwelling slabs. These lateral temperature variations in excess of 1000K may have a crucial influence on the rheology of mantle rocks. Here we use a numerical multiscale model that allows us to make predictions from first principles with no adjustable parameters on the deformation of MgO under the extreme conditions of mantle pressure, temperature and strain rate, in order to investigate the sensitivity of mantle viscosity to the temperature heterogeneities inferred from a global high resolution mantle circulation model. Our results show that under the very low strain rates of the mantle, MgO deforms mostly at low stresses (few tens of MPa) in an athermal regime, where the deformation is insensitive to both temperature and strain rate, leading to a very weak phase throughout much of the upper half of the lower mantle. In its lower half, the weak phase gives way to high material strength with thermally activated viscosities in the cold downwelling slabs, while much of the hot upwelling flow remains in the athermal regime, resulting in large lateral variations in the inferred material strength of MgO. Our results suggest the presence of large lateral viscosity variations in the deepest parts of the lower mantle, associated in particular with the graveyard of old subducted oceanic lithosphere.

Diffusion Model Coupled with the Flory–Huggins–Zuo Equation of State and Yen–Mullins Model Accounts for Large Viscosity and Asphaltene Variations in a Reservoir Undergoing Active Biodegradation

A very large monotonic variation of asphaltenes and viscosity has been measured by downhole fluid analysis (DFA) in crude oils in five-stacked sandstone reservoirs in the northern part of the Barmer Basin, northwest India, undergoing active biodegradation; each of the five-layered sand bodies shows overlaying fluid gradients with depth providing replicate validation of the measurements. Fluid data from four wells across the field shows that the gradients are uniform across the formation. The crude oil in the upper half of the oil column exhibits an equilibrium distribution of asphaltenes matching predictions of the Flory–Huggins–Zuo Equation of state (FHZ EoS) with the gravity term only using asphaltene nanoaggregates of the Yen–Mullins model. However, the bottom half of the reservoir reveals a large asphaltene gradient approximately three times larger than the equilibrium predictions from the FHZ EoS. This increase in asphaltenes creates a very large (8×) viscosity gradient and is a major production conc...

Effect of GIA models with 3D composite mantle viscosity on GRACE mass balance estimates for Antarctica

Seismic data indicate that there are large viscosity variations in the mantle beneath Antarctica. Consideration of such variations would affect predictions of models of Glacial Isostatic Adjustment (GIA), which are used to correct satellite measurements of ice mass change. However, most GIA models used for that purpose have assumed the mantle to be uniformly stratified in terms of viscosity. The goal of this study is to estimate the effect of lateral variations in viscosity on Antarctic mass balance estimates derived from the Gravity Recovery and Climate Experiment (GRACE) data. To this end, recently-developed global GIA models based on lateral variations in mantle temperature are tuned to fit constraints in the northern hemisphere and then compared to GPS-derived uplift rates in Antarctica.We find that these models can provide a better fit to GPS uplift rates in Antarctica than existing GIA models with a radially-varying (1D) rheology. When 3D viscosity models in combination with specific ice loading histories are used to correct GRACE measurements, mass loss in Antarctica is smaller than previously found for the same ice loading histories and their preferred 1D viscosity profiles. The variation in mass balance estimates arising from using different plausible realizations of 3D viscosity amounts to 20 Gt/yr for the ICE-5G ice model and 16 Gt/yr for the W12a ice model; these values are larger than the GRACE measurement error, but smaller than the variation arising from unknown ice history. While there exist 1D Earth models that can reproduce the total mass balance estimates derived using 3D Earth models, the spatial pattern of gravity rates can be significantly affected by 3D viscosity in a way that cannot be reproduced by GIA models with 1D viscosity. As an example, models with 1D viscosity always predict maximum gravity rates in the Ross Sea for the ICE-5G ice model, however, for one of the three preferred 3D models the maximum (for the same ice model) is found near the Weddell Sea. This demonstrates that 3D variations in viscosity affect the sensitivity of present-day uplift and gravity rates to changes in the timing of the ice history. In particular, low viscosities (<1019 Pas) found in West Antarctica make the mantle very sensitive to recent changes in ice thickness.

Crustal rheology control on earthquake activity across the eastern margin of the Tibetan Plateau: Insights from numerical modelling

The devastating 2008Mw7.9 Wenchuan and 2013 Mw6.6 Lushan earthquakes occurred in the Longmen Shan (LMS) fault zone, the eastern margin of the Tibetan Plateau. Seismology investigation reveals that earthquakes of Ms>5.0 are frequently recorded in the eastern Tibetan Plateau, however rare earthquake of Ms>5.0 and only a few earthquakes of 3.0<Ms<5.0 have been recorded in the Sichuan Basin. This study investigates the relation between the crustal rheology and the earthquake activity across the eastern margin of the Tibetan Plateau through 3D numerical experiments assuming visco-elasticity rheology. We setup several finite element lithospheric models based on different rheological structure. The interseismic stress accumulation processes in the models are simulated with boundary conditions of steady compressional deformation rate constrained by GPS observations. The results show that the crustal stress accumulation process across the eastern margin of the Tibetan Plateau is significantly controlled by the regional crustal rheology, which presents large horizontal viscosity variation in the lower crusts between the Tibetan Plateau and the Sichuan Basin. The stress accumulation rate in the upper crust of the eastern Tibetan Plateau is much higher than the Sichuan Basin. This stress spatial distribution explains the seismology observation that earthquakes are densely recorded in the eastern Tibetan Plateau, whereas rarely in the Sichuan Basin. The stress concentration nearly at the bottom of the entire LMS fault zone is thought to be responsible for the generation of the 2008Mw7.9 Wenchuan and 2013Mw6.6 Lushan earthquakes. There is a high possibility of future earthquakes in the seismic gaps at the southwest segment of the LMS fault zone, where it was stress accumulated but did not rupture during the last two events of the 2008 Wenchuan and 2013 Lushan earthquakes.