Subgrid Flux Research Articles

A Stable Implementation of a Data‐Driven Scale‐Aware Mesoscale Parameterization

AbstractOcean mesoscale eddies are often poorly represented in climate models, and therefore, their effects on the large scale circulation must be parameterized. Traditional parameterizations, which represent the bulk effect of the unresolved eddies, can be improved with new subgrid models learned directly from data. Zanna and Bolton (2020), https://doi.org/10.1029/2020gl088376 (ZB20) applied an equation‐discovery algorithm to reveal an interpretable expression parameterizing the subgrid momentum fluxes by mesoscale eddies through the components of the velocity‐gradient tensor. In this work, we implement the ZB20 parameterization into the primitive‐equation GFDL MOM6 ocean model and test it in two idealized configurations with significantly different dynamical regimes and topography. The original parameterization was found to generate excessive numerical noise near the grid scale. We propose two filtering approaches to avoid the numerical issues and additionally enhance the strength of large‐scale energy backscatter. The filtered ZB20 parameterizations led to improved climatological mean state and energy distributions, compared to the current state‐of‐the‐art energy backscatter parameterizations. The filtered ZB20 parameterizations are scale‐aware and, consequently, can be used with a single value of the non‐dimensional scaling coefficient for a range of resolutions. The successful application of the filtered ZB20 parameterizations to parameterize mesoscale eddies in two idealized configurations offers a promising opportunity to reduce long‐standing biases in global ocean simulations in future studies.

Interpretable Multiscale Machine Learning‐Based Parameterizations of Convection for ICON

AbstractMachine learning (ML)‐based parameterizations have been developed for Earth System Models (ESMs) with the goal to better represent subgrid‐scale processes or to accelerate computations. ML‐based parameterizations within hybrid ESMs have successfully learned subgrid‐scale processes from short high‐resolution simulations. However, most studies used a particular ML method to parameterize the subgrid tendencies or fluxes originating from the compound effect of various small‐scale processes (e.g., radiation, convection, gravity waves) in mostly idealized settings or from superparameterizations. Here, we use a filtering technique to explicitly separate convection from these processes in simulations with the Icosahedral Non‐hydrostatic modeling framework (ICON) in a realistic setting and benchmark various ML algorithms against each other offline. We discover that an unablated U‐Net, while showing the best offline performance, learns reverse causal relations between convective precipitation and subgrid fluxes. While we were able to connect the learned relations of the U‐Net to physical processes this was not possible for the non‐deep learning‐based Gradient Boosted Trees. The ML algorithms are then coupled online to the host ICON model. Our best online performing model, an ablated U‐Net excluding precipitating tracer species, indicates higher agreement for simulated precipitation extremes and mean with the high‐resolution simulation compared to the traditional scheme. However, a smoothing bias is introduced both in water vapor path and mean precipitation. Online, the ablated U‐Net significantly improves stability compared to the non‐ablated U‐Net and runs stable for the full simulation period of 180 days. Our results hint to the potential to significantly reduce systematic errors with hybrid ESMs.

Generalized eddy‐diffusivity mass‐flux formulation for the parametrization of atmospheric convection and turbulence

AbstractThe atmospheric convection is a phenomenon that has a length scale smaller than the resolution at which the state‐of‐the‐art general circulation numerical models are running, and thus the convection needs to be parametrized. In this work, a theoretical formulation for the parametrization of subgrid‐scale convection based on subgrid averaging that represents a generalized eddy‐diffusivity mass‐flux formulation is presented. The subgrid fluxes are derived by considering a decomposition of subgrid variables into convective and turbulent variables, and by assuming that the convection is modeled by round convective plumes with generic radial profiles. The main difference between our formulation and the mass‐flux formulations is that the condition of a very small fractional area occupied by the convection is replaced with the condition that the convective vertical velocity goes to the mean state far from the updraft region of the plume. The plume model can also be generalized by considering that the convective elements are energy‐consistent plumes, governed by the conservation of mass, momentum, kinetic energy, and buoyancy, which provides a physical‐based closure for the entrainment. Therefore, the closing problem of our formulation consists of the specification of the convective radial profiles and the boundary conditions at the initial vertical level. Furthermore, our formulation allows one to consider more realistic profiles for the convective variables, making the formulation suitable for the parametrization of atmospheric convection at the convective gray zone. This aspect is discussed in the last part of this work, where it is showed how the formulation can be implemented at any resolution. Moreover, in the present framework, no distinct assumptions are required for each convective type, thus facilitating the parametrization of convection in a unified way.

Subgrid Parameterizations of Ocean Mesoscale Eddies Based on Germano Decomposition

AbstractOcean models at intermediate resolution (1/4°), which partially resolve mesoscale eddies, can be seen as Large eddy simulations of the primitive equations, in which the effect of unresolved eddies must be parameterized. In this work, we propose new subgrid models that are consistent with the physics of two‐dimensional flows. We analyze subgrid fluxes in barotropic decaying turbulence using Germano (1986, https://doi.org/10.1063/1.865568) decomposition. We show that Leonard and Cross stresses are responsible for the enstrophy dissipation, while the Reynolds stress is responsible for additional kinetic energy (KE) backscatter. We utilize these findings to propose a new model, consisting of three parts, that is compared to a baseline dynamic Smagorinsky model. The three‐component model accurately simulates the spectral transfer of energy and enstrophy and improves the representation of KE spectrum, resolved KE and enstrophy decay in a posteriori experiments. The backscattering component of the new model (Reynolds stress) is implemented both in quasi‐geostrophic and primitive equation ocean models and improves statistical characteristics, such as the vertical profile of eddy KE, meridional overturning circulation and cascades of kinetic and potential energy.

Confronting the Convective Gray Zone in the Global Configuration of the Met Office Unified Model

AbstractIn atmospheric models with kilometer‐scale grids the resolution approaches the scale of convection. As a consequence the most energetic eddies in the atmosphere are partially resolved and partially unresolved. The modeling challenge to represent convection partially explicitly and partially as a subgrid process is called the convective gray zone problem. The gray zone issue has previously been discussed in the context of regional models, but the evolution in regional models is constrained by the lateral boundary conditions. Here we explore the convective gray zone starting from a defined global configuration of the Met Office Unified Model using initialized forecasts and comparing different model formulations to observations. The focus is on convection and turbulence, but some aspects of the model dynamics are also considered. The global model is run at nominal 5 km resolution and thus contributions from both resolved and subgrid turbulent and convective fluxes are non‐negligible. The main conclusion is that in the present assessment, the configurations which include scale‐aware turbulence and a carefully reduced and simplified mass‐flux convection scheme outperform both the configuration with fully parameterized convection as well as a configuration in which the subgrid convection parameterization is switched off completely. The results are more conclusive with regard to convective organization and tropical variability than extratropical predictability. The present study thus endorses the strategy to further develop scale‐aware physics schemes and to pursue an operational implementation of the global 5 km‐resolution model to be used alongside other ensemble forecasts to allow researchers and forecasters to further assess these simulations.

A Posteriori Learning for Quasi‐Geostrophic Turbulence Parametrization

AbstractThe use of machine learning to build subgrid parametrizations for climate models is receiving growing attention. State‐of‐the‐art strategies address the problem as a supervised learning task and optimize algorithms that predict subgrid fluxes based on information from coarse resolution models. In practice, training data are generated from higher resolution numerical simulations transformed in order to mimic coarse resolution simulations. By essence, these strategies optimize subgrid parametrizations to meet so‐called a priori criteria. But the actual purpose of a subgrid parametrization is to obtain good performance in terms of a posteriori metrics which imply computing entire model trajectories. In this paper, we focus on the representation of energy backscatter in two‐dimensional quasi‐geostrophic turbulence and compare parametrizations obtained with different learning strategies at fixed computational complexity. We show that strategies based on a priori criteria yield parametrizations that tend to be unstable in direct simulations and describe how subgrid parametrizations can alternatively be trained end‐to‐end in order to meet a posteriori criteria. We illustrate that end‐to‐end learning strategies yield parametrizations that outperform known empirical and data‐driven schemes in terms of performance, stability, and ability to apply to different flow configurations. These results support the relevance of differentiable programming paradigms for climate models in the future.

Numerical proof of shell model turbulence closure

The development of turbulence closure models, parametrizing the influence of small non-resolved scales on the dynamics of large resolved ones, is an outstanding theoretical challenge with vast applicative relevance. We present a closure, based on deep recurrent neural networks, that quantitatively reproduces, within statistical errors, Eulerian and Lagrangian structure functions and the intermittent statistics of the energy cascade, including those of subgrid fluxes. To achieve high-order statistical accuracy, and thus a stringent statistical test, we employ shell models of turbulence. Our results encourage the development of similar approaches for 3D Navier-Stokes turbulence.

Dynamic tensorial eddy viscosity model: Effects of compressibility and of complex geometry

A previous paper by Cimarelli et al. [“General formalism for a reduced description and modelling of momentum and energy transfer in turbulence,” J. Fluid Mech. 866, 865–896 (2019)] has shown that every decomposition of turbulent stresses is naturally approximated by a general form of tensorial eddy viscosity based on velocity increments. The generality of the formalism is such that it can also be used to give a reduced description of subgrid scalar fluxes. In the same work, this peculiar property of turbulent stresses and fluxes has been dynamically exploited to produce tensorial eddy viscosity models based on the second-order inertial properties of the grid element. The basic idea is that the anisotropic structure of the computational element directly impacts, although implicitly, the large resolved and small unresolved scale decomposition. In the present work, this new class of turbulence models is extended to compressible turbulence. A posteriori analysis of flow solutions in a compressible turbulent channel shows very promising results. The quality of the modeling approach is further assessed by addressing complex flow geometries, where the use of unstructured grids is demanded as in real world problems. Also in this case, a posteriori analysis of flow solutions in a periodic hill turbulent flow shows very good behavior. Overall, the generality of the formalism is found to allow for an accurate description of subgrid quantities in compressible conditions and in complex flows, independent of the discretization technique. Hence, we believe that the present class of turbulence closures is very promising for the applications typical of industry and geophysics.

Resolved and subgrid-scale crossing trajectory effects in Eulerian large eddy simulations of size-dependent droplet transport

We study the dispersion characteristics of slightly buoyant droplets in a turbulent jet using large eddy simulations (LES). The droplet number density fields are represented using an Eulerian approach, with the dispersed phase modelled using the Fast-Eulerian method (Ferry & Balachandar, Intl J. Multiphase Flow, vol. 27, issue 7, pp. 1199–1226, 2001) that includes the droplet rise velocity. Radial concentration profiles and turbulent concentration fluxes for droplets of different sizes are analysed to quantify the ‘trajectory crossing effect’, when relative motions between particles and turbulent eddies tend to reduce turbulent diffusion. For finer LES grid resolutions, the model captures the differential, size-based dispersion characteristics of the droplets with the transverse dispersion of the larger droplet sizes suppressed, since trajectory crossing effects are explicitly resolved in LES. We examine a similarity solution model for the size-dependent radial concentration profiles based on a modified Schmidt number derived from the theory of turbulent diffusion of particles in the atmosphere proposed by Csanady (J. Atmos. Sci., vol. 20, pp. 201–208, 1963). The results are validated with the high resolution LES data and show good agreement. Then the size-dependent Schmidt number model is reformulated as a model for unresolved subgrid-scale trajectory crossing effects and used to calculate the subgrid concentration flux in a coarse LES of a turbulent jet, with slightly buoyant droplets injected at the centreline in the self-similar region of the jet. The results are compared to a simulation with higher grid resolution and a coarse simulation with a constant Schmidt number subgrid-scale model. We find that the subgrid model enhances the prediction accuracy of the concentration profiles and turbulent concentration flux for the coarse LES.

Subgrid-scale parametrization of unresolved scales in forced Burgers equation using generative adversarial networks (GAN)

Stochastic subgrid-scale parametrizations aim to incorporate effects of unresolved processes in an effective model by sampling from a distribution usually described in terms of resolved modes. This is an active research area in climate, weather and ocean science where processes evolved in a wide range of spatial and temporal scales. In this study, we evaluate the performance of conditional generative adversarial network (GAN) in parametrizing subgrid-scale effects in a finite-difference discretization of stochastically forced Burgers equation. We define resolved modes as local spatial averages and deviations from these averages are the unresolved degrees of freedom. We train a Wasserstein GAN (WGAN) conditioned on the resolved variables to learn the distribution of subgrid flux tendencies for resolved modes and, thus, represent the effect of unresolved scales. The resulting WGAN is then used in an effective model to reproduce the statistical features of resolved modes. We demonstrate that various stationary statistical quantities such as spectrum, moments, autocorrelation, etc. are well approximated by this effective model.

Use of Neural Networks for Stable, Accurate and Physically Consistent Parameterization of Subgrid Atmospheric Processes With Good Performance at Reduced Precision

AbstractA promising approach to improve climate‐model simulations is to replace traditional subgrid parameterizations based on simplified physical models by machine learning algorithms that are data‐driven. However, neural networks (NNs) often lead to instabilities and climate drift when coupled to an atmospheric model. Here, we learn an NN parameterization from a high‐resolution atmospheric simulation in an idealized domain by accurately calculating subgrid terms through coarse graining. The NN parameterization has a structure that ensures physical constraints are respected, such as by predicting subgrid fluxes instead of tendencies. The NN parameterization leads to stable simulations that replicate the climate of the high‐resolution simulation with similar accuracy to a successful random‐forest parameterization while needing far less memory. We find that the simulations are stable for different horizontal resolutions and a variety of NN architectures, and that an NN with substantially reduced numerical precision could decrease computational costs without affecting the quality of simulations.

Observations and Simulations of Meteorological Conditions over Arctic Thick Sea Ice in Late Winter during the Transarktika 2019 Expedition

The parameterization of ocean/sea-ice/atmosphere interaction processes is a challenge for regional climate models (RCMs) of the Arctic, particularly for wintertime conditions, when small fractions of thin ice or open water cause strong modifications of the boundary layer. Thus, the treatment of sea ice and sub-grid flux parameterizations in RCMs is of crucial importance. However, verification data sets over sea ice for wintertime conditions are rare. In the present paper, data of the ship-based experiment Transarktika 2019 during the end of the Arctic winter for thick one-year ice conditions are presented. The data are used for the verification of the regional climate model COSMO-CLM (CCLM). In addition, Moderate Resolution Imaging Spectroradiometer (MODIS) data are used for the comparison of ice surface temperature (IST) simulations of the CCLM sea ice model. CCLM is used in a forecast mode (nested in ERA5) for the Norwegian and Barents Seas with 5 km resolution and is run with different configurations of the sea ice model and sub-grid flux parameterizations. The use of a new set of parameterizations yields improved results for the comparisons with in-situ data. Comparisons with MODIS IST allow for a verification over large areas and show also a good performance of CCLM. The comparison with twice-daily radiosonde ascents during Transarktika 2019, hourly microwave water vapor measurements of first 5 km in the atmosphere and hourly temperature profiler data show a very good representation of the temperature, humidity and wind structure of the whole troposphere for CCLM.

Gray Zone Partitioning Functions and Parameterization of Turbulence Fluxes in the Convective Atmospheric Boundary Layer

AbstractHere, we present the first attempt to fully represent three‐dimensional turbulence fluxes in the “Terra Incognita” or the gray zone in other words. In order to derive partitioning functions, representing the partitioning between subgrid and total fluxes, we make use of high‐resolution large‐eddy simulations (LES), which are performed with the Weather Research and Forecasting (WRF) model. LES computations are performed for various levels of convective instability, ranging from pure buoyant to strongly sheared convection. Then, the resulting reference‐LES fields are successively coarse grained from its original microscale grid spacing ( m) up to typical mesoscale grid spacings ( km). The given process is applied by means of an advanced filter, that is, the Butterworth filter. It enables a clear scale‐specific filtering that results in a more controlled energy transition from lower to higher wavenumbers, unlike the drawbacks of current filters in use. Finally, we parameterize the subgrid scale (SGS) partitioning functions of 10 SGS turbulence quantities: momentum fluxes (τij, six terms), heat fluxes (qj, three terms), and turbulence kinetic energy (k). Turbulence partitioning relations are parameterized in a scale‐aware, stability‐dependent, and height‐dependent form, using the sigmoidal Gompertz function. Thus, the new gray zone model provides a framework that bridges the mesoscale and microscale limits and that is suitable for the development of next generation three‐dimensional, multiscale turbulence parameterization methods or planetary boundary layer schemes.

Large eddy simulation of biomass gasification in a bubbling fluidized bed based on the multiphase particle-in-cell method

Numerical simulations of biomass gasification in a bubbling fluidized bed have been performed using a hybrid approach, where filtered gas-phase governing equations are solved in a Eulerian framework, and solid particles are modeled based on the multiphase-particle-in-cell method following Lagrangian stochastic parcel approach. A solver has been developed in OpenFOAM after incorporating both homogeneous and heterogeneous gasification reactions. Filtered governing equations of gas-phase variables are solved within the LES framework, where unresolved subgrid stresses and fluxes are obtained using the subgrid kinetic energy equation closure model. The solver could demonstrate transient flow characteristics of gas-solids flows in a fluidized bed, including formation, coalescence, and the eruption of bubbles at the bed surface. For different operating conditions such as steam/biomass ratio, equivalence ratio, and gasifier temperature, predicted concentrations of product gases at the gasifier outlet yield a good agreement with the experimental measurements.

Statistical Behaviors of Turbulent Scalar Fluxes in High-Pressure Turbulent Premixed Combustion in the Context of Large Eddy Simulations

ABSTRACT The effects of pressure on the physical behavior of sub-grid scalar fluxes and their modeling have been analyzed based on a-priori analysis of a Direct Numerical Simulation database of turbulent premixed Bunsen flames featuring different pressure levels. Due to an increasing ratio of the hydrodynamic to critical length scale for the occurrence of Darrieus-Landau instability, the flames become increasingly hydrodynamically unstable with increasing pressure. An increasing extent of counter-gradient transport with increasing pressure has been found, consistent with recent observations reported in the existing literature. The performance of a selection of well-known Large Eddy Simulation models for sub-grid turbulent scalar flux from literature has been assessed for a range of filter widths for different thermodynamic pressures. While some models show high correlations with the sub-grid scalar flux extracted from explicitly filtered Direct Numerical Simulation data, and a reasonable quantitative prediction of the sub-grid turbulent scalar flux magnitude for small values of filter width to flame thickness ratio, several models considerably underpredict the sub-grid scalar flux magnitudes for large filter widths. This effect becomes more pronounced for higher pressures as the length scale separation between the filter width and flame thickness increases with increasing pressure.

Analysis of the Closures of Sub-grid Scale Variance of Reaction Progress Variable for Turbulent Bunsen Burner Flames at Different Pressure Levels

The statistical behaviour and modelling of the sub-grid variance of reaction progress variable have been analysed based on a priori analysis of direct numerical simulation (DNS) data of turbulent premixed Bunsen burner flames at different pressure levels. An algebraic expression for sub-grid variance, which can be derived based on a presumed bi-modal sub-grid distribution of reaction progress variable with impulses at unburned reactants and fully burned products, has been found to be inadequate for the purpose of prediction of sub-grid variance even for the flames in the wrinkled flamelets/corrugated flamelets regime. Moreover, an algebraic model, which is often used for modelling sub-grid variance of passive scalars, has been found to significantly overpredict the sub-grid variance of reaction progress variable for all the cases considered here. The modelling of the unclosed terms of the sub-grid variance transport equation has been analysed in detail. Suitable model expressions have been identified for the sub-grid flux of variance, reaction rate contribution and scalar dissipation rate based on a priori analysis of DNS data. It has been found that the alternation of pressure does not have any significant impact on the closures of sub-grid flux of variance but a model parameter for the Favre-filtered scalar dissipation rate needs to be modified to account for the variation of pressure.

Sub-grid Reaction Progress Variable Variance Closure in Turbulent Premixed Flames

The statistical behaviour and modelling of the sub-grid variance of reaction progress variable have been analysed based on a priori analysis of Direct Numerical Simulation (DNS) data of freely propagating statistically planar turbulent premixed flames with different turbulence intensities. It has been found that an algebraic expression, which can be derived based on a presumed bi-modal sub-grid distribution of reaction progress variable with impulses at unburned reactants and fully burned products, is inadequate for the purpose of prediction of sub-grid variance. An algebraic model, which is often used for modelling sub-grid variance of passive scalars, has been found to significantly overpredict the sub-grid variance of reaction progress variable for the cases considered here. The statistical behaviours of the terms of the sub-grid variance have been analysed in detail and explained based on scaling arguments. It has been found that reaction rate and molecular dissipation contributions play leading order roles in the transport of sub-grid reaction progress variable variance and they remain in approximate equilibrium for large filter widths. Suitable model expressions have been identified for the sub-grid flux of variance, reaction rate contribution and scalar dissipation rate based on a priori analysis of DNS data.

The coupling of deep convection with the resolved flow via the divergence of mass flux in the IFS

The resolution of the European Centre for Medium‐range Weather Forecast (ECMWF) integrated forecast system (IFS) is expected to reach 5 km in the coming decade. Assumptions in the parametrization of deep convection, such as that all of the compensating environmental flow occurs in the grid column, i.e. the convective and environmental mass fluxes cancel each other in term of mass transport, have to be challenged. In this paper, we further develop the original concept of separating the convective updraught from the subsiding branch of the overturning convective circulation and apply it to the global hydrostatic equations of the IFS. In practice, this constitutes a revised convection–dynamics coupling where the mass flux subsidence of the dynamical variables is not computed locally by the convection scheme, but instead is recomputed from the revised continuity equation and is effective through the semi‐Lagrangian advection of the dynamical core. Therefore horizontal divergence/convergence is also generated in the dynamics at the top/bottom of the convective columns, thus adding to the representation of deep convection a three‐dimensional character which is not present in traditional schemes. The proposed physics–dynamics coupling is intended to be applicable to any mass flux convection scheme and within any regional or global model. We first demonstrate the accuracy of the revised physics–dynamics coupling in terms of global temperature and moisture budgets. The potential impact of the coupling on the convective organization is demonstrated for an idealized squall line case at high horizontal resolution using the small planet testbed. Model reforecasts at 9 km and 5 km resolution confirm the viability of the method in terms of forecast skill and model climate. However, the model impacts are limited as the main factor that still determines the convective stabilization and organization is the current conceptual model of subgrid mass flux which, actually, remains unchanged.

Turbulent scalar fluxes in H2-air premixed flames at low and high Karlovitz numbers

The statistical behaviour and closure of sub-grid scalar fluxes in the context of turbulent premixed combustion have been assessed based on an a priori analysis of a detailed chemistry Direct Numerical Simulation (DNS) database consisting of three hydrogen-air flames spanning the corrugated flamelets (CF), thin reaction zones (TRZ) and broken reaction zones (BRZ) regimes of premixed turbulent combustion. The sub-grid scalar fluxes have been extracted by explicit filtering of DNS data. It has been found that the conventional gradient hypothesis model is not appropriate for the closure of sub-grid scalar flux for any scalar in the context of a multispecies system. However, the predictions of the conventional gradient hypothesis exhibit a greater level of qualitative agreement with DNS data for the flame representing the BRZ regime. The aforementioned behaviour has been analysed in terms of the properties of the invariants of the anisotropy tensor in the Lumley triangle. The flames in the CF and TRZ regimes are characterised by a pronounced two-dimensional anisotropy due to strong heat release whereas a three-dimensional and more isotropic behaviour is observed for the flame located in the BRZ regime. Two sub-grid scalar flux models which are capable of predicting counter-gradient transport have been considered for a priori DNS assessment of multispecies systems and have been analysed in terms of both qualitative and quantitative agreements. By combining the latter two sub-grid scalar flux closures, a new modelling strategy is suggested where one model is responsible for properly predicting the conditional mean accurately and the other model is responsible for the correlations between model and unclosed term. Detailed physical explanations for the observed behaviour and an assessment of existing modelling assumptions have been provided. Finally, the classical Bray–Moss–Libby theory for the scalar flux closure has been extended to address multispecies transport in the context of large eddy simulations.

Application of the Eulerian subgrid Probability Density Function method in the Large Eddy Simulation of a partially premixed swirl flame

ABSTRACTA gas turbine model combustor is studied using Large Eddy Simulation with a transported Probability Density Function approach solved by the Eulerian stochastic field method. The chemistry is represented by a reduced methane mechanism containing 15 steps and 19 species while the subgrid scale stresses and scalar fluxes are modelled, respectively, via a dynamic Smagorinsky model and a gradient diffusion approximation. The test case comprises a partially premixed swirl flame in a complex geometry. Four stochastic fields are utilised in the simulations, which are performed for two different combustor operating conditions involving a stable and an unstable flame. Good agreement between the simulation and measurement data is shown in a comparison of mean velocity, temperature and species mass fraction profiles, as well as scatter plots of the instantaneous thermochemical properties. In conclusion, the predictive capabilities of the employed Large Eddy Simulation method are successfully demonstrated in this work.