Study on separation factor models for well anti-collision analysis

Anti-collision analysis has been becoming even more important in the past few years with the increasing amount of wells drilled in highly congested fields. The separation factor (SF) is a critical safety parameter to avoid wellbore collisions, where the pedal curve method (PCM) is the commonly used way for the SF calculation. However, such method may be conservative in many situations and possibly triggers an unnecessary early stop drilling due to its conservation. In this study, our work aims to evaluate and compare the different types of the SF models. In addition, an alternative and supplement method of the PCM, called as elliptic cylinder of uncertainty (ECoU) method, is proposed. The radius of the ECoU considers the distance of the intersecting points of the center-to-center line and the cylinder surfaces of both reference and offset wells. By introducing the cylinder surfaces into the SF calculation, it could extend the collision avoidance along the whole well path to reduce multiple-points collision possibility. The findings in the study help us for better understanding of the SF analysis and also show that using the ECoU indication has a great potential in fields applications as a viable method to keep the drilling operations’ safety.

Directional uncertainty with ellipsoid radius for offset well

Background
The enhancement of the oil and gas recovery in mature fields is leading worldwide companies to increase a number of wellbores in the same area. This situation makes the probability of hitting another wells to grow, potentially resulting in financial and life losses. Improving the accuracy of determining wellbore positions and keeping the safe distance among wells become important and demanding, that need more research studies on such topic. For example, in congested fields, a collision between two wells is a severe problem. To save costs and keep safety, it is crucial to increase the importance of extending and improving the production by increasing the number of available slots without the expense of significant additions to infrastructure (Poedjono 2009). For those reasons, an accurate anti-collision analysis shall be performed. The anti-collision analysis identifies the position of a planned well in relation to other nearby wells. It determines if a wellbore can be drilled safely without colliding with others. The anti-collision plan starts by determining the accurate positions of all wellbores in the area and ends with plans for well paths of future wells (Rocha et al. 2011). A typical anticollision analysis process starts by collecting the necessary data (location and position uncertainties) from offset wells in the vicinity where a well is planned to be constructed. It is also necessary to implement the position error model to estimate the planned wellbore's position uncertainties. The end product of the analysis is to determine how safe the current planned well path is. If the planned well path is not enough safe according to the requirements, the trajectory is altered until it is safe to be executed. The process is repeated for each offset well drilled in the same area.
Many companies have their own anti-collision analysis process, making the industrial standardization very difficult. In addition, mathematical and computational difficulties due to the nature of the 3D problem of the anti-collision analysis make it harder to apply available solutions. The thorny issue requires a simple representation that all drilling process actors can understand and use (Sawaryn et al. 2018). The Industry Steering Committee on Wellbore Survey Accuracy (ISCWSA), also known as SPE Wellbore Positioning Technical Section (WPTS), was created to address not only the use and accuracy of the directional drilling instruments, but also the issues related to anti-collision, well interception, and industry education, etc. They have been working through the years to promote the standardization of the anti-collision analysis.
The separation factor (SF) indicates how close two well paths are to each other. At the current state, there is no industry standard formulas to calculate the SF. However, the commonly used way to calculate a ratio between the center to center (C-C) distance of two wells and consider the wellbore position uncertainties for each center point, where the C-C distance is the one between a point in the reference well and the closest point in the offset well (Jamieson et al. 2007). The position uncertainty of a well path at a specific depth is caused due to many different error sources, and it can be geometrically expressed as an Ellipsoid of Uncertainty (EoU). The ellipsoid radius in the direction of the closest point of the offset well is commonly calculated by the pedal curve method (PCM) (Dabyah et al. 2016). The PCM calculates the size of the position uncertainty in the direction from a point in the reference well path to the other point in the offset well such that the results are dependent on the direction. Although the PCM is considered the safest method for the SF calculation, it is very conservative. In certain situations, the SF calculated by the pedal curve radius (Lawrence 1972) can erroneously indicate that two wells (reference and offset wells) have collided.

Objectives and novelty
To have better understanding of the different calculations on the SF, in our study, the three versions of the SF models 1 3 will be compared and evaluated. The SF equation provided by the ISCWSA (ISCWSA 2019) is the recommended model due to its robustness and stability. Besides it, a realistic and field-tested SF equation, provided by Equinor (a Norwegian state-owned energy company) is considered in this work. Lastly, a simple version of the equation recently mentioned by Mansouri et al. (2020) is also used. From the analysis on the different SF models, the model suggested by the ISCWSA shows the more robust and adequate in real operations, see the discussions in Sect. 6.3.
To calculate the SF, most of the existing methods consider the bit error ellipsoids of the reference and offset wells, instead of an envelop of a series of ellipsoids along the wellbore paths. It can neglect the possibility that the collisions might happen in multiple-points/sections along well paths. In our study, a new method is proposed to compensate such limitation quite well by introducing a 3D surface of cylinder uncertainties of the reference well and offset well to the SF calculation, where the collision can be avoided along the whole well path to make the collision possibility minimum.
The EoU (see more detailed descriptions in Sect. 2) at a particular depth can be derived from a 2D cut of the ellipsoid, with the center at a point in the well path and being orthogonal to the trajectory direction. A series of ellipses of uncertainty separated by a small distance in the well path, and connected in a chain, form a surface of uncertainty around the well path. In our work, this surface is called the elliptic cylinder of uncertainty (ECoU). The size of the ECoU is calculated by finding the radius of the intersecting points of the C-C line and the cylinder surfaces for both reference and offset wells. More details of this method are given in Sect. 4. Moreover, in Sect. 6, the SF values based on the ECoU are compared with the ones using the EoU. From the comparisons, it shows that in some scenarios, the ECoU can precisely indicate when the two well paths have in factor collided, being less pessimistic than the PCM. This paper first presents the theory on the EoU in Sect. 2. Then, the different SF models are summarized in Sect. 3. Section 4 presents a method to calculate the ECoU. After it, Sect. 5 lists the main procedure for the SF calculation. Following it, Sect. 6 shows the main results in our study to compare the SFs by using the different SF models and with the use of EoU/ECoU.

Ellipsoid of uncertainty
For a given survey tool, the multiple factors may lead to bit position errors. They are evaluated at the particular survey station and used to determine position uncertainties. Errors from the survey stations accumulate along a wellpath as measure depth (MD) increases, and the contribution of all individual errors is summed to determine the total uncertainty of the wellbore position. These uncertainties are expressed as a covariance matrix ( COV NEV ) which describes the error along each of the coordinate axes, North-East-Vertical (NEV) coordinate. The computational steps and mathematical details of calculating COV NEV can be found in (ISCWSA 2019). Its mathematical expression is given as below: To calculate the SF, the radius of the bit ellipsoid needs to be considered. An uncertainty ellipsoid centered at the point c can be expressed below, given in Friendly et al. (2011): In the literature, there are two main methods to determine the radius of the EoU, called as pedal radius and ellipsoid radius, respectively.
The pedal radius is defined as an orthogonal intersection between the C-C line and a line tangent to the ellipsoid, see Fig. 1 (double arrow in purple). It can be also referred to as the ellipsoid's projected dimension onto the C-C line (the line between two crosses). The method to calculate the pedal curve radius was presented by Bang et al. (2020) to calculate the root squared of the multiplication between the covariance matrix and the unit tangent vector of ⃗ r . The P 1 (r) and P 2 (r) are the individual directional uncertainties of the reference and the offset well, respectively, as pedal radius: The ellipsoid radius is a radius of the ellipsoid along the C-C line or the distance from the point of the well path to the point of the ellipsoid surface where the C-C line intersects it, see Fig. 1 (double arrow in yellow). Equation (4) to determine the directional uncertainties E 1 (r) and E 2 (r) as ellipsoid radius was presented in Bang et al. (2020): Figure 1 shows the difference between the ellipsoid radius and the pedal radius, where the ellipsoid radius has the higher opportunity to underestimate the uncertainty for some relative orientations of the ellipsoids. Because of it, the ellipsoid radius is not currently implemented in most commercial anti-collision applications. Currently the pedal radius is the default method for most well-planning software (ISCWSA 2013). (1)

Separation factor
The SF is a concept widely used by the Petroleum Engineering industry to determine if the distance between two wells is safe. Unfortunately, despite many efforts, there is no a standard way to determine it. Different companies are using different versions, creating confusions to influence the safety of the drilling operations and possibly leading to catastrophic results. Although, as mentioned, different SF models are used, the ratio factor is the type of the equation recommended by the ISCWSA, since it represents a close relationship between the safe distance and the probability of the well crossing (Bang et al. 2020). The ratio-type is commonly a ratio between the C-C distance and the directional uncertainties. In the following, three different SF models are presented that are often used in research work or companies for anticollision analysis.
The first SF model is to calculate the SF from the general ratio-type SF equation. It has been used extensively by many companies and planning software. Recently it was used by Mansouri et al. (2020) for the well path optimization design. Here SF s represents the SF calculated from such model: However, the SF s is the least stable of the SF models. At shallow depths, the calculated SF might be very large due to the small directional uncertainty D . When the wells are near, Eq. (5) might be misleading and potentially dangerous, since it does not include the dimensions (well radius) (ISCWSA 2013). In order to improve it, for instance, Equinor, a Norwegian state-owned energy company, has used a slightly different version of the SF equation. The primary SF ratio is modified to incorporate the radius of both reference and offset wells. The corresponding SF equation becomes: Moreover, the SF model proposed by the ISCWSA in Sawaryn et al. (2019), is widely considered as the robust one, given below: where the safety margin S m defines the minimum acceptable separation during the wells' design while the project-ahead uncertainty pa is partially correlated with the projection from the survey depth to a distance beyond the bit. The recommended values by Sawaryn et al. (2019) are S m = 0.3 m and pa = 0.5 m and they have much more impact on shallow depths where the significance of the parameters is increased. Equation (7) is with an effort to standardize the anti-collision analysis across the different companies. The WPTS believes that the standardization will clarify the expectations and requirements and support consistency between operators, service providers, and regulators and increase planning and operational efficiencies and reduce the burden of training (Sawaryn et al. 2019). Table 1 summarizes these three SF models. For the first SF model ( SF s ), it uses a linear summation of the term 1 and 2 to calculate D . As pointed out in Bang et al. (2020), the linear summation method to combine the individual uncertainty produces pessimistic (smaller) SF results. For SF e and SF wpts calculations, the directional uncertainty D is calculated by combining the individual covariance matrices of the points in the reference and offset well. For all of three SF models, it would be recommended to use the pedal radius ( j = P j , j = 1, 2 ) for the directional uncertainty calculation, see Table 1. A collision between wellbores, as already mentioned, can be determined by the SF value. If the value equals or is less than SF critic , a technical collision most likely happens. In Bang et al. (2020), the authors recommend SF critic = 1 to be used, although many companies have adopted through Fig. 1 Example of the ellipsoid radius and the pedal curve radius the years a stop drilling criterion of SF critic = 1.25 , or even SF critic = 1.5 for the most conservative ones.

Elliptic cylinder of uncertainty
In this section, a new approach to calculate the uncertainty radius j , j = 1, 2 based on the ECoU is presented. It can be further used for the SF calculation given in Sect. 3. A series of ellipsoids of uncertainty along the well path with the center at the survey stations can form a surface that envelops the ellipsoids, see Fig. 2 (blue envelop for offset well, gray envelop for reference well). The resultant surface is an elliptic cylinder with principal radii size and posture defined by the rotation and the ellipsoids' principal radii size. In Liu (2019), the mathematical equation of such cylinder is defined, that is expressed in the local coordinating system [x, y, z], given as below: The details of the cylinder surface calculation is given in Appendix A.
The surface of the ECoU is very intricate, where the size of its principal radii and the posture are changing with the measured depth. To solve it numerically, one point defined as ⃗ p c , where the C-C line touches the cylinder surface is considered, given as below: The point that the C-C vector line intersects the surface of the cylinder ⃗ p c is found when the unit directional vector of ⃗ v equals to the unit directional vector of C-C vector line ⃗ r .  Fig. 2 Example of two plotted cylinders for reference and offset well Therefore, the cylinder surface radius for both reference and offset wells can be calculated, respectively, by the equations below: Figure 3 shows a comparison among the pedal radius distance, the ellipsoid radius distance, and the cylinder surface distance. More discussions about their main differences will be given in Sect. 6.2. Once C 1 (r) and C 2 (r) are calculated, the SF will be determined with respect to the different SF models given in Table 1 by setting j = C j .

Implementation routine
An algorithm to calculate the SF is summarized below, from getting the survey measurements to outputting the anti-collision report listing the separation factors per measured depth. The steps are as follows:

Survey Measurements Conversion
The global coordinates are calculated from the MD, Inclination, and Azimuth survey measurements. The method implemented uses a vector equation to calculate the North, East, and Vertical positions presented by Sawaryn and Thorogood (2005). 2. Bit Error Model As described by ISCWSA (2017a), the purpose of the bit error model is to evaluate the effects of the various physical factors which lead to errors in the survey measurements. Those physical factors need to be identified for a given surveying tool (i.e., MWD tools), and effects on each of the measurements (MD, Inclination, and Azimuth) must be determined. The combination of all individual errors gives the uncertainty in the position. 3. Position Uncertainty After the error sources with its magnitudes, weighting function, and propagation type are defined, the size of the uncertainty error in the NEV coordinate system is then calculated. The core error formula was presented in ISCWSA (2017a), see the more detailed descriptions in Appendix B. 4. Covariance Matrix Covariance matrix ( COV NEV,j ) for one particular error source (j) can be easily calculated for random error sources or systematic and global error sources in ISCWSA (2017a). Per error source, a covariance matrix (3 × 3) is generated in the NEV coordinate for each survey station. The summation of all co-variance matrices generated that way will produce the final co-variance matrix COV NEV per survey station. 5. Ellipsoid of Uncertainty Since the covariance matrix is a positive semi-definite matrix and assuming the positional error has a normal distribution, it can be represented by a uncertainty ellipsoid, as shown in Eq.
(2) to calculate the EoU. 6. Directional Uncertainty The equations for calculating the pedal radius, the ellipsoid radius and the ECoU radius are given in Eqs.
(3), (4) and Eqs. (10)-(11), respectively. 7. Separation Factor Three equations for the SF calculation are presented in Table 1. The one could select either P j , E j or C j to get j .

Reference and offset wells
The wellbores were selected from the ISCWSA clearance table, which was designed to help users implement and verify their anti-collision algorithms. The clearance table in ISCWSA (2017b) has 1 reference well and 11 offset wells. For this work, only four offset wells were selected. Table 2 shows the list of wellbores. The configuration data used in this work are given in Table 3. Figs. 4 and 5 give a good overview of the reference well and offset wells.

C-C, pedal, ellipsoid, and cylinder distances
The C-C distance is calculated using the minimum distance scanning method (Rocha et al. 2011). Figures 6 and 7 comparing the four types of distances. The pink distance is the C-C distance between points in the reference well and the offset well. The distance in blue is the one between cylinder surfaces. Finally, the green and the yellow lines are the distance between ellipsoids when calculating 1 and 2 by the pedal radius method and the ellipsoid radius method, respectively.
The yellow curve is the distance between ellipsoids calculated using the ellipsoid radius method. It uses Eq. (4) to determine the point where the C-C line "touches" the ellipsoid boundary. This curve is plotted to contrast with the cylinder surface distance, which also uses the concept of the C-C line intersecting the surface. The cylinder surface can be treated as the sum of ellipses of uncertainty throughout the whole well path. As displayed in Figs. 6 and 7, the distance between the survey point and the cylinder surface considers all the ellipsoid of uncertainties through the whole wellbore. While both the pedal curve radius and the ellipsoid radius methods are only considering the ellipsoid associated with the actual survey point. In extreme cases, using only the pedal curve radius or the ellipsoid radius method, can make the SF more optimistic than in reality.
As it can be seen in Figs. 6 and 7, when the wellbores are almost parallel to each other, as the North100 offset well is parallel to the reference well, the distances are almost the same, except for the ellipsoid radius method (yellow curve) that after 1100 m (MD) starts to present the separation from the other distances. For the offset wells: East10 and East20, it is observed negative distance. It is possibly due to the sum of individual directional uncertainties 1 and 2 that are more significant than the C-C distance. This situation possibly indicates that either the ellipsoids or the cylinder surfaces from both the reference and offset wells are "touching" each other.

3
The pedal curve radius method's limitation was pointed out in Jamieson et al. (2007). It depends on the ellipsoids' posture and how close they are. The PCM can be very pessimistic (low SF values) as if the surfaces were touching each other. As shown in Fig. 1, the ellipsoids are not touching each other, which could still be considered a safe situation (depending on the confidence level chosen for each ellipsoid). However, the SF would be potentially less than 1, unnecessarily alerting for a collision situation that did not happen. On the other hand, when using the cylinder surface radius, it shows precisely that, although the wellbores are reasonably close to each other, the uncertainty areas are not touching each other, and it is not considered a technical collision. The SF using the cylinder surface distance would be a small value, maybe below most of the operator's safety Considering the situation presented in Figure 1, both the SF calculated by cylinder surface distances or the ellipsoid radius distance method, would be potentially precisely the same. It does not mean that solutions are similar. In Fig. 3, the SF calculated by the cylinder surface distance and the ellipsoid radius distance would be dramatically different depending on the direction of the C-C line. The cylinder's SF would be much smaller than the one produced by the ellipsoid radius distance or even the pedal radius distance.

Separation factors
The ISCWSA Eq. (7) is considered throughout this work as the base equation and is effectively being used to compare and verify the results using the cylinder surface distances and the results using the pedal radius distance. Equation (6) is brought to this discussion as an example of a SF used daily in the field in real operations. At last, the version of the SF ratio: SF s as an alternative way of calculating the importance factor, was utilized in this work to compare its performance facing the other two. 1 All three basic equations have different ways of calculating the D variable. In order to compare the results, a common scaling factor, k = 3.058 , equivalent to a confidence level of 97.5% , is applied to the following simulations.

Comparisons: SF s , SF e and SF wpts
The results by Eq. (5) are more conservative (smaller value) than the others. Also, Eq. (5) does not consider the radius of the reference and the offset wellbores, which makes it less robust and safer than the other two.
As can be seen from Figs. 8 and 9, the differences between Eqs. (6) and (7) for the pedal radius distance are marginal while Eq. (5) presents in general a smaller safety factor than the other two equations.

Comparisons: D selection
In this case, Eq. (7) is used for the SF calculation. For the D selection, three options are considered here for comparisons. The first one is calculated from Eq. (3) by using pedal radius (red one); the second one is calculated from Eq. (4) with the use of ellipsoid radius (blue one); the third one is calculated from Eqs. (10)-(11) based on the cylinder surface radius (orange one).
As it can be seen in Figs. 10 and 11, the cylinder surface case was matching almost perfectly to its respective equation until the geometry and the attitude of the wellbores (reference and offset) towards each other started to play a significant role in the calculation. Except in the offset well North100, which is parallel to the reference well, the offset wells' safety factors started to present differences after a certain depth. In the offset well East10, the curves until around 2000 m (MD) were matching very closely. After that depth, the distances, as observed in Fig. 6, also started to reduce. The pedal radius distance reduces faster than the cylinder surface distance. The ellipsoid radius method took longer to identify the risk than the other two reflecting in the separation factors at the same depth.

Advantages and limitations
As displayed in Fig. 1, the pedal radius method can be very pessimistic, presenting a situation where the sum of the 1 and 2 are greater than D, while the cylinder would present a smaller value for both terms. For some situations, especially when the C-C line is almost orthogonal to both reference and offset well paths, the cylinder surface distance and the ellipsoid radius will be the same. However, when the angle between the C-C line and the path is less than 90 • , the cylinder will present a much more accurate safety factor than the other two. Figure 3 shows an example of the advantages of using the cylinder distances instead of the pedal or ellipsoid radius methods. The findings show that using the cylinder distances have a good potential as an alternative measure factor for anti-collision analysis.
When it has analyzed the relationship between the reference and the Angular offset well, it was observed that when calculating the distances using the pedal radius, the distance between ellipsoids is reduced to less than 0 while the cylinder surface distances are still in the positive side. For certain situations, both the ellipsoid radius method and the cylinder surface distance method share the same weakness: it can be very optimistic, leading to the possible collisions. However, the ellipsoid radius method results are quite different from the cylinder surface results in the Angular Offset well case. The SF calculations based on cylinder surface, pedal radius or ellipsoid radius have different limitations and disadvantages. It would recommend to use two or three radius for the SF calculations. It would be important to attract more attentions on such topic to introduce more accurate ways for the SF calculations.

Summary and conclusions
From the SF models analysis, it is clear that the model provided by the ISCWSA is the most robust and less prone to numerical errors than the other two. It demonstrates in this work that the SF s model is more pessimistic, displaying much smaller SFs than the others when using the PCM. It also does not use both wellbores' radius in the calculation, making it dangerous to be used. The biggest difference between the Equinor and ISCWSA's SF equations is the additions of the safety margin and the projected-ahead uncertainties introduced by the ISCWSA. Those terms consider a minimum separation between the wellheads and the projection uncertainties, making the equation from the ISCWSA more robust and adequate to be used in real operations.
By comparing the three methods to calculate the individual directional uncertainties, 1 and 2 , it is clear that the pedal curve radius is the most conservative, displaying the smallest SFs as the two wells, the reference, and the offset, are getting close to each other. As already explained, the reason is due to the very nature of how the pedal curve is derived, and it can cause an earlier stop drilling. The cylinder surface method is presented as a valid alternative to the pedal curve, although it faces some issues similar to the ellipsoid radius method. To overcome this issue, the well planning engineer would have to determine the smallest distance between the two ellipsoids to ensure the position uncertainties do not overlap each other. With more wellbores being drilled in the same area in the past few years, it is now more than important to define an industry-standard way to defined how safe is to drill in the area, avoiding costly well collisions, without forgetting its precision. The elliptic cylinder of uncertainty can be considered as a good alternative to the pedal curve. The main conclusions are summarized below: • The SF model provided by the ISCWSA is the most robust and less prone to numerical errors; • The pedal curve radius is very conservative for certain scenarios; • The elliptic cylinder of uncertainty can be considered as a good alternative to the pedal curve.

Appendix A
The cylinder is expressed in the local coordinating system [x, y, z], given as below: The transformation matrix B can be found from the relation between the eigenvectors of COV NEV and the transformation matrix T: The transformation matrix T is: As pointed by Liu (2019), the projection curve equation of the tangent curve between the elliptic cylinder surface and the ellipsoid on the plane xy is: where F is a (2 × 2) matrix with its terms defined as follows: The term D i,j are defined as follows: and C j : If the matrix G = F −1 , as mentioned by Liu (2019) the radius of the principal axes and the rotation angle of the cross section ellipse of the error elliptic cylinder, similar to the ellipses of uncertainty principal axes equations, are: The cylinder surface was found using the basic elliptic cylinder parametric equations: where ∈ [0, 2 ], v ∈ [0, ΔMD] for a very small ΔMD. The elliptic cylinder surface is found for the [x, y, z] coordinate system. To change the coordinating system to nev coordinate system, the transformation matrix T needs to be multiplied to [x, y, z] cylinder coordinates. The new parametric equations of the etcou are as follows: D i,j = B i,j − B i,3 C j , for i = 1, 2, 3 j = 1, 2 R cyl,1 = √ G 1,1 cos 2 + G 1,2 sin 2 + G 2,2 sin 2 R cyl,2 = √ G 1,1 sin 2 − G 1,2 sin 2 + G 2,2 cos 2 = 1 2 tan −1 2G 1,2 G 1,1 − G 2,2 .

Appendix B
The core equation used for the error evaluation is: where represents the error source (magnetometer, accelerometer etc.), i is the index for the particular error source considered and i is the magnitude of the uncertainty for the ith error. The weighting function p is a 3x1 vector with respect to the survey measurements for depth, inclination, and azimuth, or In addition, d d is the effect of the survey errors in measured depth, inclination and azimuth on the wellbore position in NEV-axes. This is expressed as an 3x3 matrix ISCWSA (2017a). Wellbore positions can be calculated using numerous methods. The d d matrix depends on the interval between two survey stations, where dr can be expressed as Δr k for the displacement between survey station k − 1 and k, and Δr k+1 for the displacement between stations k and k + 1 . Equation (19) can be written with respect to Δr k and Δr k+1 .
In the definition of the ISCWSA model ISCWSA (2017a), the balanced tangential method is used to calculate the dr k dp k . It provides the equations below for displacement between two survey stations k − 1 and k in the NEV-axes.
The final 3x3 matrix equations of dΔr k dp k and dΔr k+1 dp k can be expressed be as follows: (21) e i,k = i dΔr k dp k + dΔr k+1 dp k p k i .
sin I k−1 cos A k−1 + sin I k cos A k sin I k−1 sin A k−1 + sin I k sin A k cos I k−1 + cos I k Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creat iveco mmons. org/ licen ses/ by/4. 0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Funding The authors received no specific funding for this work.
Availability of data and material The authors will provide the data and material used in the study through the open-access software.