Comparison of an impec and a semi-implicit formulation for compositional reservoir simulation

PROCESS SYSTEMS ENGINEERING

Comparison of an impec and a semi-implicit formulation for compositional reservoir simulation

B. R. B. Fernandes^I; A. Varavei^II; F. MarcondesIII, ^*; K. Sepehrnoori^IV

^ILaboratory of Computational Fluid Dynamics, Federal University of Ceará, Brazil

^IICenter for Petroleum and Geosystems Engineering, The University of Texas at Austin, USA

^IIIDepartment of Metallurgy and Materials Science and Engineering, Federal University of Ceará, Brazil. E-mail: [e-mail]

^IVCenter for Petroleum and Geosystems Engineering, The University of Texas at Austin, USA

ABSTRACT

In compositional reservoir simulation, a set of non-linear partial differential equations must be solved. In this work, two numerical formulations are compared. The first formulation is based on an implicit pressure and explicit composition (IMPEC) procedure, and the second formulation uses an implicit pressure and implicit saturation (IMPSAT). The main goal of this work is to compare the formulations in terms of computational times for solving 2D and 3D compositional reservoir simulation case studies. In the comparison, both UDS (Upwind difference scheme) and third order TVD schemes were used. The computational results for the aforementioned formulations and the two interpolation functions are presented for several case studies involving homogeneous and heterogeneous reservoirs. Based on our comparison of IMPEC and IMPSAT formulations using several case studies presented in this work, the IMPSAT formulation was faster than the IMPEC formulation.

Keywords: Compositional reservoir simulation; Segregated formulation; IMPEC; IMPSAT; Finite-volume method.

INTRODUCTION

Several formulations have been developed for solving the governing partial differential equations arising from modeling fluid flow for compositional simulations in porous media. In general, formulations can be classified as Implicit Pressure Explicit Composition (IMPEC), Implicit Pressure and Saturation (IMPSAT), or the Fully Implicit Method (FIM). Phase saturation is defined as a volume of phase per void space available for fluid flow. The IMPEC formulation has the lowest cost in terms of computational time per time-step. However, due to the higher degree of explicitness in the calculation of composition, this formulation cannot use large time-steps when compared to the FIM and IMPSAT approaches. The IMPSAT method can handle larger time-steps, compared to the IMPEC approach. It is also less expensive in terms of computational time, per time-step, than the FIM approach. Additionally, the IMPSAT approach is more stable than the IMPEC formulation due to the reduction in the degree of explicitness. Also, according to Cao (2002), it has good performance compared to other FIM approaches, since saturations are much more coupled than compositions. Therefore, although the saturation calculation involves solution of a linear system of equations for the IMPSAT formulation, the larger time-steps used by the formulation compensate the overall cost and render a CPU time reduction compared to IMPEC approaches.

In this work, an IMPSAT formulation proposed by Watts (1986) was implemented into the UTCOMP simulator. UTCOMP was developed at the Center for Petroleum and Geosystems Engineering at The University of Texas at Austin for the simulation of enhanced recovery processes. The UTCOMP simulator is a multiphase/multi-component compositional equation-of-state simulator, which can handle the simulation of several enhanced oil recovery processes. The original numerical procedure of the UTCOMP simulator is an IMPEC formulation based on Àcs et al. (1985).

Several procedures for solving pressure and saturation implicitly (IMPSAT) have been proposed in the literature (Branco and Rodrigues, 1996; Kendall et al., 1983; Spillette et al., 1973). However, the IMPSAT approach proposed by Watts (1986) and adopted in this work has similar features to the original IMPEC formulation of the UTCOMP simulator. For instance, the IMPSAT formulation solves simple sets of linear systems for both pressure and saturations, while only one flash procedure is performed per time-step, thus allowing the calculations per time-step to be much faster than the other IMPSAT and IMPEC approaches that need iterations, thus performing flash and solving the conservation equations until convergence. The approaches that do not iterate in a time level are called one-iteration formulations; therefore, this work is based on a comparison of two one-iteration approaches: an IMPEC and an IMPSAT formulation. Although the IMPEC and IMPSAT formulations implemented and used in this work are not new, to the best of our knowledge this is the first time that these formulations are compared for three phase hydrocarbon flow simulations (oil, gas, and a second liquid hydrocarbon phase). The second liquid hydrocarbon phase is important for CO₂ injection processes, where the CO₂ and some light components tend to form a CO₂-rich phase. Another important feature included in this work is the use of a high-resolution TVD scheme to approximate the fluxes for the Watts' formulation. Also, only few performance results are shown in the literature for this formulation. Some of these results can be found in Haukås (2006). In this work, results are compared in terms of volumetric production rates, saturation fields, and CPU time. For most investigated cases, the formulation was able to use larger time-steps than the IMPEC formulation, achieving the same results.

PHYSICAL MODEL

The Watts' formulation is basically an adaptation of the method of Spillette et al. (1973) combined with the formulation of Ács et al. (1985), where the volume error constraint is added to pressure and saturation equations in order to use just only one flash calculation per time-step. In this section, we show the molar balance equations, the pressure equation, and the new saturation equations that are included in the original formulation of Ács et al. (1985).

If advection is the only transport mechanism involved, the molar balance equations according to Chang (1990) are given by

where N_k is the moles of component k, V_b is the bulk volume, x_kj is the mole fraction of component k in phase j, ξ_j is the molar density, respectively, q_k is the molar rate of component k through the well, N_c is the number of hydrocarbon components, N_c+1 denotes the water component, k_rj and µ_j are the relative permeability and viscosity of phase j, respectively, is the absolute permeability tensor, and Φ_j is the hydraulic potential of phase j, which is defined by

where P is the pressure of the oil phase, ρ_j is the mass density of phase j, g is the gravity, D is the depth, which is positive in the downward direction, and P_crj is the capillary pressure between phases j and r.

Fluid phase equilibrium between the hydrocarbon phases is considered (water is not considered in any flash calculation). This assumption considers that the chemical potential of all phases are the same. This can be expressed in terms of the equality of the fugacities (f) of the phases, which can be stated as follows:

Fugacity and PVT properties (molar density, compressibility factors, and volume derivatives) are evaluated in this work using the Peng-Robinson Equation of State (EOS) (Peng and Robinson, 1978). The flash procedure used considers a fixed and known pressure, temperature, and global compositions (isothermal flash) in order to evaluate the phase compositions and fluid properties. Further details of this procedure can be found in Perschke (1988). A volume-shift approach based on the work of Jhaveri and Youngren (1988) is also available for liquid density correction.

Two phase stability test algorithms are implemented in the UTCOMP simulator: the stationary point location method (Michelsen, 1982) and the Gibbs free energy minimization algorithm, which is similar to the Trangenstein (1987) method and was modified by Perschke (1988) to deal with three hydrocarbon phase equilibrium. In general, as commented by Perschke (1988), the stationary method is faster than the Gibbs free energy method; therefore, the stationary method was used in this work. The flash calculation used in UTCOMP is a combination of the Accelerated Successive Substitution (ACSS) method (Mehra et al., 1983) with the modified version of the Gibbs free energy minimization method (Perschke, 1988). At the beginning of the flash procedure, we use the ACSS method in order to provide a reasonable initial estimation, and then we switch to the Gibbs free energy minimization method in order to accelerate the convergence. The switching criterion to change from one method to another is given by Chang (1990) as:

where, the superscript 2 denotes the oil phase. The switching criterion (ε_swi) equal to 0.01 as suggested by Chang (1990) is used.

The pressure equation used for both formulations is based on the volume constraint proposed by Ács et al. (1985). The pressure equation is obtained from the equality between the formation pore volume (V_p) and the total fluid volume (V_t):

where the pore volume is given by:

where φ_o is the porosity at the reference pressure (P_f), and c_f is the rock compressibility.

Taking the derivative of Eq. (5) with respect to time, applying the chain rule to the right-hand side, substituting Eqs. (1) and (6), and dividing all by the bulk volume, we obtain

where N denotes derivative evaluated by holding the number of moles constant.

The saturation equation is solved implicitly only for the IMPSAT formulation. The approach used here is described by Watts (1986). By definition, the saturation of phase ℓ (S_ℓ), as stated before, is the ratio of the phase volume to the volume of the pore:

where V_ℓ is the volume of phase ℓ, which is a function of pressure and number of moles. Equation (8) can be written as:

Taking the derivative of Eq. (9) with respect to time, applying the chain rule on the right-hand side and substituting Eq. (1) and dividing it by the bulk volume, we obtain:

Comparing Eqs. (7) and (10), we can observe that saturations and pressure equations have similar formats. As suggested by Watts (1986) and Spillette et al. (1973), the phase velocity that is needed in Eq. (10) should be evaluated as a function of total velocity. According to these authors, if this approach is taken, the segregated solution of pressure and saturations will be equivalent to a fully implicit procedure in terms of pressure and saturations. Next, we present the expression for the total velocity using the idea given in Watts (1986), but now considering a full permeability tensor. Extending the steps of Kaasschieter (1999) for a three-phase flow, or even four-phase flow, yields:

where

and,

where, is the total fluid velocity, λ_m is the m-th phase mobility which is defined as a ratio of relative permeability and viscosity of the m-th phase, and f_j is called the fractional flow, which is defined as above.

Substituting Eq. (11) into (10) yields:

Equation (15) is the final saturation equation in terms of total velocity. The next section is devoted to showing the numerical discretization applied to the above equation.

APPROXIMATE EQUATIONS

In order to obtain an approximate equation for the saturation of phase ℓ, we will integrate Eq. (15) over the control volume of Figure 1 and time.

The final form of the saturation equation for the three-dimensional control volume as shown in Figure 1, is given by:

where is the total volume derivative with respect to total number of moles of component k, and is a semi-implicit velocity at the east interface, which is given by:

where,

The semi-implicit velocity, at the other interfaces, is given by equations similar to Eq. (17). Newton's method is used to treat the non-linearities involved in Eqs. (16) and (17) on account of the implicit evaluation of the relative permeabilities and capillary pressures.

Performing a similar procedure for the saturation equation, we obtain the approximate mole balance equations, in terms of the semi-implicit velocity, as given below:

The final form of the pressure equation is obtained (considering only the x direction for simplicity) as:

Further details regarding the approximate equations for the pressure and moles can be obtained in Chang (1990).

Fig. 2 shows the flow chart of the procedure used in the UTCOMP simulator for the implementation of the Watts' formulation. As we can see in this figure, pressure is evaluated first; then, an iterative procedure is used to evaluate the saturations.

The mobilities, the densities, and the mole fraction of each component in each phase at each gridblock interface are evaluated using a one-point upwind interpolation function. A total variation diminishing (TVD) method using the Koren's flux limiter (Koren, 1993; Liu et al., 1995; Fernandes et al., 2013) was also implemented for solving the saturation equations in the relative permeabilities and phase compositions. In this situation, as a non-linear term is included in the relative permeability due to the flux-limiter, it was treated semi-implicitly using the idea presented by Rubin and Blunt (1991), in which the higher order terms are only considered in the independent term of the linear system. For the one-point upwind implementation, the mobilities are evaluated as:

Solving the problem using the one-iteration IMPEC approach is straightforward since only one set of linear equations (for pressure) needs to be solved and the total number of moles is computed semi-implicitly to conserve the moles. A flash calculation is performed to obtain phase compositions and phase mole fractions. Finally, the water and the hydrocarbon saturations, respectively, are evaluated as:

and

where L is the phase mole fraction obtained in the flash in the absence of water.

In this work, relative permeability is modelled by the Corey's model (Corey, 1986) and capillary pressure is modelled according to Chang (1990) as:

where C_pc and E_pc are user-input parameters, σ is the interfacial tension between the two phases and denotes a normalized saturation.

For both approaches used in this work, one needs to solve the same number of variables. All variables are listed in Table 1. Also, Table 1 shows some functional relations for some of the variables used in this work.

The phase appearance for the IMPEC and IMPSAT approaches developed in this work is treated only by the phase stability test, as described earlier. Also, the phase disappearance for the IMPEC approach is based on the phase stability test. However, for the IMPSAT approach, because phase saturation can be computed as zero or negative during the solution of Eq. (16), we have to couple the stability test with another approach. In this case, if the solution of Eq. (16) produces a negative saturation, we set that saturation to zero. When a negative saturation is calculated, it means that the phase disappears for that time-step. Setting saturation to zero is a common approach used for some fully implicit approaches; see Coats (1980), for instance.

RESULTS AND DISCUSSION

In this section, we investigate the numerical solutions, as well as the performance in terms of CPU time of the Watts' formulation implemented in this work and the original IMPEC approach of the UTCOMP simulator. The comparison studies were carried out by empirically setting the maximum allowable time-step for various case studies, which did not produce oscillatory results for both investigated formulations.

The first case investigated is CO₂ injection in an isotropic heterogeneous reservoir. All reservoir data used for this case are shown in Table 2. The water is not considered in flash calculations and is not injected in the reservoir for this case. The water mole numbers are estimated by density and initial saturation.

The components, the initial fluid compositions, and the injected fluid compositions are shown in Table 3.

The absolute permeability fields in the x and y directions are shown in Fig. 3. In order to better visualize the variation of the permeability fields, Figs. 3b and 3c show two different zooms of the whole scale presented in Fig. 3a.

We compare the results of the Watts' formulation with the original IMPEC formulation of the UTCOMP simulator. The results in terms of oil and gas production are presented in Figs. 4a and 4b, respectively. From Fig. 4, we can observe a good agreement between the original IMPEC formulation and the Watts' formulation implemented.

In this case, water saturation (volumetric fraction) is lower than the so-called connate saturation (0.25) and hence water cannot flow through the rock pores. Initially, only oil and a mixture of CO₂ and light hydrocarbons injected (gas) are present in the reservoir. As a result, CO₂ displaces oil first and forms a third non-aqueous phase. This phase is referred to as a second liquid phase (not considering the aqueous phase) (see Figure 5). It is important to mention that this behavior is only possible by considering three-phase flash calculations. Also, the present case showed a slow convergence when we tried to run using a two-phase flash (only oil and gas). This reinforces the importance of considering the three-hydrocarbon flash calculations. The results for the second liquid front at 1844 days, obtained with the IMPEC and the Watts' formulation, are compared in Figs. 5. As can be seen in this figure, a good agreement at the two fronts is verified.

The time-steps used by the original UTCOMP approach (IMPEC) and the Watts' approach are shown in Fig. 6.

As we can observe in Fig. 6, the time-steps used by the Watts' formulation were larger than those of the IMPEC formulation for the entire simulation. Due to the computational time spent for the solution of the linear system equations for the saturations, the Watts' formulation is more expensive per time-step. However, the large time-steps used by the Watts' formulation allowed this formulation to be less expensive in terms of CPU time than the IMPEC formulation. This fact can be seen in Table 4, which shows the total CPU time used by both formulations. From this table, we observe that Watts' formulation is about two times faster than the IMPEC formulation. We also check the formulations' performances by doubling the number of gridblocks in each direction. As we can see in Table 4, the speed-up ratio of the Watts' formulation is improved when the number of gridblocks is doubled in each direction.

Case 2 is similar to Case 1, but we replaced the upwind interpolation function by a third-order TVD interpolation function. Figs. 7a and 7b show the volumetric rate of oil and gas, respectively, and Fig. 8 presents the second liquid saturation front at 1844 days. Once again, the Watts' formulation results are in good agreement with the IMPEC formulation. Despite the large difference when comparing the production rates for TVD solution (Figure 7) with that of upwind solution (Figure 4), the phase behavior obtained was similar for the two cases (the third phase does appear).

The time-steps used in this case study for the IMPEC and the Watts' formulations are shown in Fig. 9. Although the time-step pattern presented in Fig. 9 is different from the ones shown in Fig. 6, using the upwind scheme, approximately the same speed-up ratio of the Watts' formulation compared to the IMPEC formulation for Case 1 was obtained as shown in Table 5.

Case 3 refers to a WAG (water-alternating gas) process in a heterogeneous reservoir. For this case study, in each cycle, we first inject CO₂, and then water. Each fluid is injected at a fixed pressure over the course of ten days. The total simulation time is 280 days, which corresponds to fourteen cycles. Tables 6 and 7 present the reservoir data set and the initial and the injected fluid compositions, respectively, employed for this case study. Two hydrocarbon phases were considered in this case and water still is not accounted for in the flash calculations. The absolute permeability field in the x and y directions is shown in Figure 10, while the permeability in the z direction is constant and equal to 7.89x10^-15 m² (8 mD). For this case study, the upwind scheme was used to obtain the solution for both formulations. Figures 11a and 11b present the oil and gas volumetric rates, respectively, and Fig. 12 shows the oil saturation front at 254 days using a 100x100x5 Cartesian grid.

The time-steps used by both formulations are shown in Fig. 13 and the total CPU time and the speed-up ratio are presented in Table 8. From the figures, we can observe that the average time-step used by Watts' formulation is approximately three times larger than that used by the IMPEC approach. This average time-step results in a speed-up ratio of 1.7, as we can see in Table 8. In general, at the beginning of each cycle, the time-step size is very close to the minimal time-step designed for the whole simulation. This approach should be in favor of the IMPEC formulation.

The fourth case study is another WAG process, but now a homogeneous reservoir is tested. Capillary pressure is now included in order to check the algorithm performance when this physical phenomenon is relevant. Except for absolute permeabilities, all of the previous data presented in Tables 6 and 7 were used. The absolute permeabilities in the x and y-directions are set to 1.97x10^-13 m², and the one in the z-direction is equal to 9.87x10^-14 m². We use an 80x80x5 grid, and the process is simulated for 100 days. The oil and gas rates are presented in Figs. 14a and 14b, respectively. Once again, we can observe a good match between the results of the two approaches compared.

The time-step comparison is shown in Fig. 15. From this figure, we can verify that the maximum time-step employed by the Watts' formulation is about two times larger than the one used by the IMPEC approach. It is important to mention that the time-step ratio, aforementioned, is relevant for only the second and third cycles, as we can verify from Fig. 15.

Table 9 presents the CPU time and the speed-up ratio obtained for case 4. From this table, we can verify that the performance of Watts' formulation is worse compared to the previous cases, but it still performs better than the IMPEC approach. Further investigation of the current Watts' implementation needs to be performed when still more complicated physical parameters are involved.

The last case study refers to a 2D gas flood with twenty-five components. Only two-hydrocarbon phases and an immobile aqueous phase are considered. The purpose of this case is to see how Watts' approach will perform with a large number of components by looking at the impact of expensive flash calculations over the simulation performance. Tables 10 and 11 show the reservoir data and the initial and the injected fluid compositions, respectively. It is worthwhile to mention that all the components equal or higher than C₂₅₊ have identical physical properties. The main goal here was to verify the performance of the compared approaches with a very large number of components.

Figure 16 shows the oil and gas volumetric rates. From this figure, we can see a good match for oil and gas volumetric rates for both approaches.

The time-steps used by both approaches are shown in Fig. 17. From this figure, we can verify that Watts' formulation was able to handle this case study with several components using much large time-steps compared to the IMPEC formulation.

Table 12 presents the total CPU time and speed-up ratio for both approaches. We also show the speed-up ratios when gridblocks in each direction are doubled. From this table, we can infer that Watts' formulation is about 2.25 times faster than the IMPEC approach.

This suggests that, for a large number of components, the Watts' formulation should be used instead of the IMPEC approach. However, when the number of gridblocks is increased, the speed-up ratio decreased to 1.93. Although the Watts' formulation performance decreased when the gridblocks are increased, this formulation is still about 2 times faster than the IMPEC approach.

CONCLUSIONS

In this work, we implemented the Watts' formulation for compositional reservoir simulation using Cartesian grids. This formulation was included into the UTCOMP simulator. Also, two interpolation functions for evaluating the physical properties at each interface of the control volume were implemented: UDS and a third-order TVD scheme. At least for the case studies investigated, the speed-up ratio of the Watts' formulation did not change when the TVD and UDS schemes were used. For most of the case studies tested, the Watts' formulation implemented in this work was around two times faster than the original IMPEC formulation of the UTCOMP simulator. We also verify, for some case studies, that the performance of Watts' formulation persists when the mesh is refined. It was also confirmed that the new implemented formulation is more efficient than the IMPEC for cases with a large number of components.

NOMENCLATURE

Greek Letters

Superscripts

Subscripts

ACKNOWLEDGMENTS

The authors would like to acknowledge the Abu Dhabi National Oil Company for the financial support for this work. We would like to thank Dr. Chowdhury K. Mamum for his comments on this manuscript. Also, the third author would like to thank the CNPq (The National Council for Scientific and Technological Development of Brazil) for financial support through the grant No. 305415/2012-3. Finally, we would like thank the ESSS (Engineering Scientific Software Simulating) for proving the Kraken^® to pre and post-processing the results.

Submitted: November 1, 2013

Revised: January 16, 2014

Accepted: February 21, 2014

• Ács, G., Doleschall, S. and Farkas, E., General purpose compositional model. SPE Journal, 25(4), p. 543-553 (1985).
• Branco, C. M. and Rodriguez, F., A Semi-implicit formulation for compositional reservoir simulation. SPE Advanced Technology Series, 4(1), p. 171-177 (1996).
• Cao, H., Development of Techniques for General Purpose Simulators. Ph.D. Thesis, Stanford University (2002).
• Chang, Y.-B., Development and Application of an Equation of State Compositional Simulator. Ph.D. Dissertation, The University of Texas at Austin, Austin Texas (1990).
• Chang, Y.-B., Pope, G. A. and Sepehrnoori, K., A higher-order finite-difference compositional simulator. Journal of Petroleum Science and Engineering, 5(1), p. 35-50 (1990).
• Coats, K. H., An equation of state compositional model. SPE Journal, 20(5), p. 363-376 (1980).
• Corey, A. T., Mathematics of immiscible fluids in porous media. Water Resources Publication, Littleton, CO (1986).
• Fernandes, B. R. B., Marcondes, F. and Sepehrnoori, K., Investigation of several interpolation functions for unstructured meshes in conjunction with compositional reservoir simulation. Numerical Heat Transfer Part A, Applications, 64(12), p. 974-993 (2013).
• Jhaveri, B. S. and Youngren, G. K., Three-parameter modification of the Peng-Robinson equation of state to improve volumetric prediction. SPE Reservoir Engineering, 3(3), p. 1033-1040 (1988).
• Kaasschieter, E. F., Solving the Buckley-Leverett equation with gravity in a heterogeneous porous medium. Computational Geosciences, 3(1), p. 23-48 (1999).
• Kendall, R. P., Morrell, G. O., Peaceman, D. W. and Watts, J. W., Development of a multiple application reservoir simulator for use on a vector computer. SPE Middle East Oil technical conference. SPE, Manama, Bahrain (1983).
• Koren, B., A Robust Upwind Discretization Method for Advection, Diffusion and Source Terms. In: C. B., Vreugdenhil, and B., Koren, (Eds.), Numerical Methods for Advection-Diffusion Problems. Notes on Numerical Fluid Mechanics, vol. 45, p. 117-138, Vieweg, Braunschweig (1993).
• Liu, L., Delshad, M., Pope, G. A. and Sepehrnoori, K., Application of higher-order flux-limited methods in compositional simulation. Transport in Porous Media, 16(1), p. 1-29 (1994).
• Lohrenz, J., Bray, B. G. and Clark, C. R., Calculating viscosities of reservoir fluids from their compositions. Journal of Petroleum Technology, 16(10), p. 1171-1176 (1964).
• Mehra, A. R., Heidemann, R. A. and Aziz, K., An accelerated successive substitution algorithm. Canadian Journal of Chemical Engineering, 61(4), p. 590-596 (1983).
• Michelsen, J. L., The isothermal flash problem. Part I. Stability, Fluid Phase Equilibria, 9(1), p. 1-19 (1982).
• Peng, D. Y. and Robinson, D. B., The characterization of the heptanes and heavier fractions for the GPA Peng-Robinson programs. Gas Processors Association (1978).
• Perschke, D. R., Equation of State Phase Behavior Modelling for Compositional Simulator. Ph.D. Thesis, The University of Texas at Austin, Austin, Texas (1988).
• Runbin, B. and Blunt, M. J., Higher-order implicit flux limiting schemes for black-oil simulation. SPE Symposium on Reservoir Simulation, SPE, Anaheim, USA, p. 219-229 (1991).
• Spillette, A. G., Hillestad, J. G. and Stone, H. L., A high-stability sequential solution approach to reservoir simulation. Fall Meeting of the Society of Petroleum Engineers of AIME. SPE, Las Vegas, Nevada, USA (1973).
• Trangenstein, J. A., Customized minimization techniques for phase equilibrium computations in reservoir simulation. Chemical Engineering Science, 42(12), p. 2847-2863 (1988).
• Watts, J. W., A compositional formulation based on pressure and saturation equations. SPE Reservoir Engineering Journal, 1(3), p. 243-252 (1986).

*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  14 Nov 2014
• Date of issue
  Dec 2014

History

• Accepted
  21 Feb 2014
• Received
  01 Nov 2013
• Reviewed
  16 Jan 2014