Research Article | Open Access
The Application of Stefan Problem in Calculating the Lateral Movement of Steam Chamber in SAGD
It is extremely important to monitor the development of steam chambers in reservoirs in the process of SAGD (Steam Assisted Gravity Drainage). By analyzing the temperature data of monitoring wells, the extending velocity and direction of steam chamber can be obtained, which helps to regulate and improve production effectively. Based on Stefan’s inverse problem, a mathematical model is established in this paper for acquiring the extending velocity of steam chamber by using the temperature data of monitoring wells. Mathematical methods are utilized to solve the model such as Fourier transformation as well as an analytic solution using the moving velocity of steam chamber’s edges. Finally, sensitivity analyses concerning some factors such as the formation heat transfer coefficient, the temperature of steam chamber, and the rate of temperature rise of monitoring wells are made. Field application showed that the model can be used to calculate the moving velocity in real time and determine the boundary scope of the steam chambers easily.
Butler [1, 2] first proposed SAGD being applied to heavy oil reservoirs. Two horizontal wells a few meters apart are drilled on the bottom of heavy oil reservoirs. Initially, steam is circulated in both injector and producer. Heavy oil around the pair of wells is heated to decrease bitumen viscosity. The viscosity of bitumen drops several orders of magnitude so that the heavy oil could flow. The period of circulating heating is quite short. After that, steam is only injected into the upper well to heat the bitumen. With the gravity and the vertical density difference between gas and liquid, heavy oil can be produced from the lower horizontal well on the bottom. The steam which is continuously injected into the formation fills the underground space left by the oil extraction, thus forming a widening steam chamber. SAGD has been proved to be one of the effective techniques of heavy oil exploitation after nearly 30 years of field application.
Steam chambers are initially formed where the steam is injected. It is a small cavity and fast lateral extension at an early stage. When steam is continuously injected, the chamber growth takes an important role. The extending velocity of steam chamber is lower than the vertical velocity of steam chamber. After reaching the interlayer on the top, steam chamber begins to extend horizontally into the formation.
Based on these factors (e.g., the operating pressure of steam chamber, the steam oil ratio, and the development characteristics of steam chamber), the startup stage and the production stage of SAGD should be subdivided into six detailed stages: the uniform pressure cycle, the uniform increasing pressure cycle, the early SAGD production, the high-pressure production, the stable production, and failed production. Sun et al.  put forward the transition conditions of each stage in SAGD. They thought the shapes of the steam chamber in different periods are important parameters in judging transition conditions of SAGD production in different stages.
In the oil field practice, the theoretical development state of the steam chamber is often predicted by the theoretical derivation and the numerical simulation. These theoretical predictions such as the theoretical derivation and the numerical simulation are mainly based on the formation characteristics, crude oil characteristics (e.g., the porosity, the relative permeability, the initial oil saturation, and the irreducible oil saturation), and the injection-production parameters (e.g., the steam injection pressure, the steam injection temperature, the steam quality, and the steam oil ratio). The process from the known conditions belongs to the direct problem solution in physical problems. Most related researches made by previous scholars about the theoretical derivation and the numerical simulation of steam chamber development belong to the direct problem.
In terms of theoretical derivation, Butler  first focused on the lateral extension of steam chamber. His theory of slope gives the calculation formula concerning the velocity and the shape of steam chamber at the interface when the steam chamber begins to extend laterally. But according to the formula, the level velocity of the interface tends to be infinite near the top of the oil layer which obviously does not match the actual data. With the additional consideration of operating pressure, Ito and Ipek  gave a formula to simulate the boundary movement of the rising steam chambers. They concluded that the higher the operating pressure, the faster the rising velocity of steam chamber. Wei et al.  gave a prediction mode in which the movement velocity of steam chamber is the relation between the reservoir properties and injection parameters. Most of the previous studies of steam chamber development are mainly based on the hypothesis of Butler, in which the boundary of steam chamber is considered as fixed and its movement velocity is considered as constant. In the actual production process of SAGD, the heavy oil at the interface is heated to reduce the viscosity. Under the gravity, the heavy oil flows into the production wells and the interface of steam chamber is continuing to push into the deep formation. The whole movement velocity is not constant. A majority of theoretical derivations about the steam chamber development are based on the hypothesis of the symmetrical steam chamber boundary, which results in a huge difference between the theoretically predicted development form of the steam chamber and the actual monitored one. Only a few studies considered the boundary of steam chamber as moving. Pooladi-Darvish et al.  calculated the shape of steam chamber and the oil production rate over time under the influence of the moving boundary in cylindrical coordinates. Their results are consistent with the experimental data. On the basis of Pooladi-Darvish et al.’s work, Mohammad  deduced a heat conduction equation in one-dimensional coordinates considering the effects of moving boundary on the steam chamber. He also discussed the influence of the Rayleigh number on the movement and the shape of steam chamber as well as the oil production rate.
For the past 30 years, with the progress and wide spread of the commercial software of numerical simulation, a large number of numerical simulations are applied to the practice of SAGD. The numerical simulation process is usually time consuming. Moreover, the simulation results often depend on the accuracy of the input parameters. But the input parameters are hard to obtain since the reservoir itself is anisotropic. These lead to unsatisfactory results. Zhu et al.  got the simulated drop curves of the temperature and the pressure during SAGD shut-in period using numerical simulation. Based on the simulated curves and the size of the development of steam chamber under different parameters, a regression equation is obtained. They compared the theoretical drop curve with that obtained from SAGD shut-in period and used the regression equation to determine the size of the development of steam chamber. Their theory is also built on the hypothesis that steam chambers develop evenly. So the diversity of the development of the steam chamber limits the application of their theory.
The actual development state of steam chamber is monitored by following some parameters such as the micro seismic observation, the formation pressure, and the formation temperature. The direction and the overall shape of steam chamber development can be analyzed by the changes of these parameters. By analyzing the changes and using the micro seismic observation, the development of steam chamber can be obtained. Results gained through micro seismic observation are accurate. But this method is not cost-effective and cannot achieve the real-time monitoring of steam chamber. In order to understand the situation of underground steam chambers and the effect of steam injection in real time, monitoring wells are employed to monitor the changes of pressure and temperature on the site. In the process of SAGD production, pressure tests of horizontal wells are made during the shut-in period. Pressure well-testing method can be used to obtain the overall size of the steam chamber, but it is unable to identify the specific shape [10, 11]. Due to the limited application of the micro seismic observation and the pressure well-testing, the temperature data from monitoring wells are often adopted to monitor real-time development of steam chambers in practice. Common tools for monitoring the steam chamber temperature include the thermocouple, the optical fiber, and the downhole electronic thermometer temperature measuring system. In the monitoring, every 05 m is set as an observation point. After a continuous monitoring for a long time, the temperature curves of monitoring wells can be drawn. Then, through the slope and the form of the curves, we can have the general knowledge of the development situation of steam chambers. The early analysis of temperature monitoring data on the site is mainly relying on experiences in combination with other geological data. If an upward or downward inflection point exists in the formation temperature curve, it usually means that the physical properties of reservoir change greatly at the inflection point. As shown in Figure 1, an inflection point appears in the formation temperature curve at the well depth of 216 m. It can be confirmed from Figure 2 (the logging curve) that there exists an interlayer. The differences of physical properties between the interlayer and the oil layer cause the dramatic change in the temperature curve (Table 1).
Due to the anisotropy of reservoir physical properties and the different density and viscosity between heated fluid and steam, boundaries of steam chambers are not uniform. In the SAGD production site, the temperature monitoring curves present a variety of shapes. In order to know the development of steam chamber more accurately, some researchers introduced the heat conduction theory to the field temperature data analysis. Carslaw and Jaeger  proposed a mathematical theory in which the interface is fixed and heat flows into the semi-infinite formation. When steam chamber moves at a constant velocity, its front temperature distribution can be acquired through Laplace mathematical change. Their results are consistent with the front temperature of steam chamber proposed by Butler. Based on the theory of Carslaw and Jaeger, Birrell and Putnam  used thermocouple data above the steam chamber to draw the inverse conjugate error function (ICEF) of dimensionless temperature versus depth and the natural logarithm of the dimensionless temperature versus depth. When the rising velocity of steam chamber changes, the inverse conjugate error function (ICEF) curve and the natural logarithm curve of dimensionless temperature will intersect at different point along the depth. The rising velocity of steam chamber can be obtained by calculating the distance between intersections. When the steam chamber stops moving upward, the inverse conjugate error function (ICEF) and the natural log of dimensionless temperature will intersect at one point along the depth. This point is the top of steam chamber. Birrell and Putnam mainly focused on using the chart to determine the rising velocity of steam chamber but rarely covered the moving velocity of the horizontal extension of steam chamber.
2. The Stefan Problem in SAGD Process
In 1889, Stefan put forward a mathematical problem to calculate the ice melting interface and the temperature change of the ice surface over time while studying melting ice. Afterwards, those problems (e.g., the melting and solidification of the polar ice at the poles, the solidification of casts, the solidifying of building components, and the freezing of food) belong to mathematical problems of solid-liquid’s moving interface location and are referred to as Stefan problem.
Stefan problem includes direct Stefan problem and inverse Stefan problem. Direct Stefan problems are divided into two categories: one is that the temperature field distribution and the heat flow at the boundary are calculated according to the change of the phase interface position and the terminal temperature of the known solid; the other is that the change of phase interface, the temperature field distribution, and the heat flow at the boundary are calculated according to the known terminal temperature of solid and Stefan conditions. Inverse Stefan problem is that the position change of the phase-change boundary is calculated according to known solid terminal temperature and the heat flow at the terminal of solid. In the actual production of SAGD, it is important to determine the location and moving velocity of interface. The process of calculating the steam chamber moving location and velocity is part of an inverse Stefan problem.
A discrete Tikhonov regularization method was adopted by Wei and Yamamoto  to determine the basic solution of one-dimensional heat conduction equation moving boundary. Li and Wei  gave a multilayered medium one-dimensional heat conduction equation by using Cauchy data to determine the basic solution of moving boundary. In addition to Tikhonov regularization method, the iterative regularization method can be used to reconstruct the unknown boundary through the initial guess of the shape of boundary and by using minimization to calculate the deviational functional between boundary data and measurement data and to determine the boundary shape. Embedding methods of different areas were adopted by Badea and Daripa  to determine the shape of interface. They considered a new interface location to replace the originally unknown interface and converted nonlinear problems to linear problems in the new expanding areas. Usually this kind of methods is very time consuming due to iteration. In order to quickly solve this problem, a finite element technique is presented to determine the phase-change interface position of two-dimensional melting problem . This solution is based on the error minimization between the measured temperature and simulated temperature of sensor. Fredman  presented a direct method which was called line method to determine the moving boundary of the inverse heat conduction problem.
In this paper, according to the temperature monitoring data of monitoring well, Stefan inverse problem is solved analytically by Fourier transform and other mathematical methods. The real-time moving velocity of steam chamber boundaries can be calculated in real time when steam chamber extends laterally.
Most existing studies of steam chamber development assume that the shape of steam chamber is an inverted triangle. The edges of steam chamber are stable. They considered that -direction is the direction normal to the steam chamber interface and the velocity is parallel to -direction. According to those assumptions, the temperatures of monitoring wells do not reveal the real situation of the development of the steam chamber at the same depth. The depth of temperature of monitoring well should be which is lower than that of the steam chamber, as shown in Figure 3. The difference of depth between monitoring well and steam chamber is reducing with lateral extension of steam chamber. In Figure 3, is the depth and is the time; is the position of steam chamber at the time under the depth of ; is the initial distance between monitoring well and horizontal well. is the angle between the steam chamber interface and horizontal line.
The result is caused by ignoring the fact that the edges of steam chamber are unstable. There exists steam fingering phenomenon on the borders of steam chamber . When steam is injected into reservoir, instability at the displacement front evolves in the form of steam fingers due to viscosity contrasts. Through numerical simulation, it is found that the scope of the steam fingering is approximately between 0.1 and 10 cm. Irani  further pointed out that the stability of the interface of steam chamber on both sides is influenced by both the temperature and the velocity of the moving interface of steam chamber.
The injected steam from reservoir and bitumen have different temperature. The variation of temperature and the instability involved in the displacement processes are referred to as thermoviscous fingering (TVF). In practice, the displacement process (steam displaces bitumen) takes place in environments with nonuniform permeability. Steam chamber extension is proportional to permeability. The velocity of interface in permeable layers has been observed to be faster and the mechanisms of fingering are more complex in heterogeneous media. The thermoviscous fingering on both sides of steam chamber is distributed horizontally and the shape of fingering is rectangular . In the micron scale, the interface of steam and bitumen is not sloping but vertical.
In order to reduce the difference of depth between the temperature of monitoring well and the development of steam chamber, we assume that -direction is horizontal direction and velocity is parallel to -direction in this paper.
There are many factors influencing the effect of SAGD implement. In order to use monitoring data to calculate the moving velocity of steam chamber and make sure of the accuracy of calculation, we make the following basic assumptions:(1)The reservoir is isotropic and direction is the same; temperature at any point within the steam chamber is the same.(2)The heat conduction is the only heat transfer mode in the front edge of the steam chamber.(3)The monitoring data is the real reflection of the change of steam chamber at the same depth.(4)The moving velocity of steam chamber interface at a certain period of time is certain.
4. Mathematical Formulation and Solution
The formation is only heated through heat conduction between monitoring well and steam interface. Heat transfer ahead of the steam chamber can be formulated as the following: where is the formation temperature, is a position at some point between the front of steam chamber and monitoring well, and is thermal diffusion coefficients.
In our model, the shape of steam chamber is assumed as arbitrary, as illustrated in Figure 4. is monitoring well and is the measured temperature distribution curve of monitoring well. The original point is at monitoring well.
The velocity of steam chamber at the interface is assumed as constant within a time period and the position of interface can be formulated as the following:where is the velocity of the advancing front of steam chamber and is the initial distance between monitoring well and steam chamber at .
In the border, there are
We introduce two parameters ( and ) and those parameters can be formulated as the following:
Because (9) is a nonhomogeneous equation, according to the principle of superposition, (9) can be decomposed into a nonhomogeneous equation and a homogeneous equation. The new homogeneous equation is formulated as the following:
The solution of (12) is as follows:
Assuming (14) is the solution of (9),where is the undetermined coefficient of the solution regarding the moving boundary. Equation (9) is translated into Fourier series formula and merged similar terms of , . The expression of is as follows:
Another new homogeneous equation of (9) is formulated as the following:
Equation (16) is a one-dimensional nonhomogeneous linear differential equation. After solving the equation, we can get the following:
According to (5), we can finally determine the temperature distribution of the front of steam chamber:
We introduce into (19); then,
In the above equation, is the heat capacity of oil, is the density of heavy oil, and is heat conductivity coefficient.
When the steam interface is pushing forward deep into the formation, steam is condensing in the interface and heat is passed to the reservoir deep in the formation through the steam-liquid interface to make the crude oil in front of steam chamber heat up to the temperature of steam chamber. According to the principle of conservation of energy, there is the following:
In the above equation, is the average velocity of the advancing front of steam chamber when the interface of steam chamber reaches the monitoring well.
5. Solving Steps
(1)The production time is divided into series periods; represents each period of time. is linear regression at each time period; is the velocity of the advancing front of steam chamber and is a constant at each time period; there is .(2)According to (21), the temperature distribution of the front edge of the steam chamber is calculated at each time period.(3)According to (24), iterative solution is done to calculate the value of at each time period.(4)The relationship between the moving rate of steam-liquid interface and time at different depth is determined and the diagram of the changes of steam chamber interface over time is driven.
6. Sensitivity Analysis
It can be seen from Figures 5 and 6 that the influence of the coefficient of thermal conductivity and the advancing velocity of steam chamber on the distribution of the front edge of steam chamber is big. The larger the coefficient of thermal conductivity, the more uniform the temperature distribution of the front steam chamber and the smaller the reducing rate of temperature; the faster steam chamber advance, the greater the temperature change of the front steam chamber; the slower the advancing velocity of steam chamber, the slower the dropping rate of the temperature of the front steam chamber. These temperature changes of front steam chamber caused by the advancing velocity of steam chamber agree with Butler’s theory. It can be seen from Figures 7, 8, and 9 that the advancing velocity is affected not only by the coefficient of thermal conductivity and the temperature of steam chamber, but also by the temperature-rising rate of monitoring wells. The effect of the coefficient of thermal conductivity on the moving velocity of steam chamber is mainly in the early stage. As time goes by, the difference of the moving velocity caused by the coefficient of thermal conductivity gradually decreases. The higher the temperature of steam chamber, the faster the moving velocity of steam chamber. Therefore, increasing the temperature of the steam chamber helps to improve the moving velocity of steam chamber. As can be seen from Figure 9, when the daily rising rate of the temperature of monitoring point is greater than 0.1 degrees, the moving velocity of steam chamber is proportional to the temperature-rising rate of monitoring wells at the same depth; when it is less than 0.1 degrees, the influence of temperature-rising rate of monitoring wells at the same depth on the moving velocity of steam chamber is minimal. The development nearly stops at later stage.
7. Instance Analysis
In Xinjiang Feng Cheng oil field in China, the viscosity of super heavy oil is big and the formation energy is low. In the production of steam stimulation, the daily capacity at the early stage is high but declines quickly. FHWX1P well group is a group of SAGD test wells in Feng Cheng oil field and the well group is preheated by steam for 3 months before production. As shown in Figure 10, around the injection-production well group, the monitoring wells FZI104, FZI107, FZI108, and FZI109 are arranged to monitor the changes of formation temperature. Because the instrument FZI108 is damaged, its temperature monitoring data cannot be used.
The temperature profile of FZI109 at the depth of 180 m mainly presented a linear distribution. The temperature profile of FZI104 at the same depth is divided into three parts and the temperature can be viewed as linear rise in each part, as shown in Figure 11. The regressions between the temperature and production time are made, as shown below.
The temperature regression of FZI104 is as follows:
The temperature regression of FZI109 is as follows: .
The temperature distribution data between the front steam chamber and monitoring wells can be calculated. Subsequently, the relation between the moving velocity of steam chamber and the change of time can be obtained through iterative calculation. According to our calculation, the lateral moving velocity of steam chamber is 0.0033~0.1867 m per day.
It can be seen from Figure 12 that, in the process of SAGD production, oil production is proportional to the moving velocity of steam chamber and the oil production cycle coincides with the fluctuation cycle of the moving velocity steam chamber. According to the calculated lateral moving velocity of steam chamber, we can finally calculate the lateral profile map of steam chamber between monitoring wells and the injection-production well group. As shown in Figures 13 and 14, the horizontal distance of steam chamber development in our calculation is 9.44 m and the steam chamber has extended 11.06 m from interpretation of microseism. The profile shape of steam chamber is similar to the temperature monitoring data in Figure 15. The good agreement between data from microseism and calculation proved that our formula has practical meaning. Therefore, through the temperature data from monitoring wells, to a certain extent, the development situation of steam chamber can be known.
(1)A new Stefan model for calculating the velocity of moving boundary is established in this paper. Based on the observed temperature data of observation wells, the real-time moving velocity of steam chamber can be calculated and the scope of the lateral development of steam chamber can be determined.(2)The value of time step directly affects the accuracy of the lateral movement of steam chamber; when dividing time periods, the periods of time in which the temperature of observation wells presents a linear relationship should be divided together as a time period.(3)When the injection pressure and injection temperature are kept constant, the moving velocity steam chamber declines with time and the declining range is related to the formation coefficient of thermal conductivity and the temperature of steam chamber; the larger the thermal conductivity and the temperature of steam chamber, the smaller the declining range. The daily temperature-rising range of monitoring wells is proportional to lateral moving velocity of steam chamber.
On the assumption that is the solution of (9), we can derive from the formula that Convert in (9) to form of Fourier series:The substitution of (A.2) and (A.3) into (9) yieldsWe merge the similar terms of in (A.4):The solution of theabove equations can be obtained as follows:
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
The authors would like to thank NSFC (National Natural Science Foundation of China) for supporting this paper through a project (no. 51304164), research on the pressure dynamics of multiple-acidized-fractured horizontal wells in fractured-vuggy carbonate formations. The paper also was financially supported by a basic research project under Grant no. 2015JY0132 from Science and Technology Department of Sichuan Province and the 973 Program (National Basic Research Program) of China under Grant no. 2014CB239205. The authors would also like to thank the reviewers and editors whose critical comments were very helpful in preparing this paper.
- R. M. Butler, G. S. McNab, and H. Y. Lo, “Theoretical studies on the gravity drainage of heavy oil during in-situ steam heating,” Canadian Journal of Chemical Engineering, vol. 59, no. 4, pp. 455–460, 1981.
- R. M. Butler and D. J. Stephen, “The gravity of steam heated to parallel horizontal wells,” Journal of Canadian Petroleum Technology, vol. 20, no. 2, pp. 90–96, 1980.
- X.-G. Sun, W.-J. He, X.-B. Hu et al., “Paramaters optimization of different production stages by dual horizontal well SAGD Process for super heavy oil Reservoir,” Xin Jang Petroleum Geology, vol. 33, no. 6, pp. 697–699, 2012.
- R. M. Butler, “Rise of interfering steam chambers,” Journal of Canadian Petroleum Technology, vol. 26, no. 3, pp. 70–75, 1987.
- Y. Ito and G. Ipek, “Steam-fingering phenomena during SAGD process,” in Proceedings of the SPE International Operations and Heavy Oil Symposium, SPE-97729-MS, Calgary, Canada, November 2005.
- S.-L. Wei, L.-S. Cheng, S.-J. Huang, and W. Huang, “Steam chamber development and production performance prediction of steam assisted gravity drainage,” in Proceedings of the SPE Heavy Oil Conference Canada, pp. 24–33, Alberta, Canada, June 2014.
- M. Pooladi-Darvish, W. S. Tortike, and S. M. Farouq Ali, “Non isothermal gravity drainage under conduction heating,” in Proceedings of the Annual Technical Meeting of CIM, PETSOC-94-65, Calgary, Canada, June 1994.
- R. H. Mohammad, Mathematical modeling of SAGD and VAPEX processes [M.S. thesis], University of Calgary, Alberta, Canada, May 2008.
- L. Zhu, F. Zeng, and G. Zhao, “An empirical correlation of steam chamber size and temperature fall off in the early-period of SAGD process,” in Proceedings of the SPE Heavy Oil Conference, SPE-165389-MS, Calgary, Canada, June 2013.
- J. Rabb and C. Palmgren, “Pressure transient analysis in SAGD,” Journal of Canadian Petroleum Technology, vol. 44, no. 9, pp. 16–22, 2005.
- A. Zamani, B. R. James, and R. H. Hite, “Use of pressure transient analysis for monitoring SAGD steam chamber development,” in Proceedings of the SPE Heavy Oil Conference, SPE-165489-MS, Calgary, Canada, June 2013.
- H. S. Carslaw and J. C. Jaeger, Conduction of Heat in Solids, Oxford Science Publications, 1959.
- G. E. Birrell and P. E. Putnam, “A study of the influence of reservoir architecture on SAGD steam chamber development at the underground test facility,” in Proceedings of the Canadian International Petroleum Conference, PETSOC-2000-104, Calgary, Canada, June 2000.
- T. Wei and M. Yamamoto, “Reconstruction of a moving boundary from Cauchy data in one-dimensional heat equation,” Inverse Problems in Science and Engineering, vol. 17, no. 4, pp. 551–567, 2009.
- Y. S. Li and T. Wei, “Identification of a moving boundary for a heat conduction problem in a multilayer medium,” Heat and Mass Transfer, vol. 46, no. 7, pp. 779–789, 2010.
- L. Badea and P. Daripa, “A domain embedding method using the optimal distributed control and a fast algorithm,” Numerical Algorithms, vol. 36, no. 2, pp. 95–112, 2004.
- Y. Ruan and N. Zabaras, “An inverse finite-element technique to determine the change of phase interface location in two-dimensional melting problems,” Communications in Applied Numerical Methods, vol. 7, no. 4, pp. 325–338, 1991.
- T. P. Fredman, “A boundary identification method for an inverse heat conduction problem with an application in ironmaking,” Heat and Mass Transfer, vol. 41, no. 2, pp. 95–103, 2004.
- D. R. Gotawala and I. D. Gates, “Flow and energy dynamics at the edges of steam-assisted gravity drainage chambers,” in Proceedings of the Canadian International Petroleum Conference, Calgary, Canada, 2008.
- M. Irani, “Understanding the steam-hammer mechanism in SAGD wells,” SPE Journal, vol. 18, no. 6, pp. 1181–1201, 2013.
- M. Sajjadi and J. Azaiez, “Thermo-viscous fingering in heterogeneous media,” in Proceedings of the SPE Heavy Oil Conference, SPE-157722-MS, Calgary, Canada, June 2012.
Copyright © 2015 Xiong Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.