**1. Introduction**

Stimulated by developments in computer power and mechanics theories, Computational Fluid Dynamic (CFD) methods were widely applied in ocean engineering to study the hydrodynamic performance of marine structures. Generally, computation domains containing target structures need to be discretized into multiple mesh elements to solve differential equations numerically. In ship maneuvering simulations, the presence of CFDbased tools enables researchers to obtain more precise hydrodynamic coefficients compared with that of traditional empirical methods, which is beneficial to establish individualized and accurate mathematical maneuverability models for different objectives. As the most widely used steering devices for ships, rudders play an important role in the maneuvering of ships [1], and the accuracy of lift coefficients has a great impact on motion predictions. Many CFD simulations were conducted to study rudder hydrodynamics (Liu et al. [2], Badoe et al. [3], Van Nguyen and Ikeda [4]), but few concentrate on mesh generations

**Citation:** Lu, S.; Liu, J.; Hekkenberg, R. Mesh Properties for RANS Simulations of Airfoil-Shaped Profiles: A Case Study of Rudder Hydrodynamics. *J. Mar. Sci. Eng.* **2021**, *9*, 1062. https://doi.org/ 10.3390/jmse9101062

Academic Editors: Carlos Guedes Soares and Serge Sutulo

Received: 31 July 2021 Accepted: 5 September 2021 Published: 28 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

for rudders. In this paper, the study on mesh properties for RANS simulations of rudder hydrodynamics is performed to present detailed insight into meshes around rudders.

A good mesh is fundamental to the convergence and the accuracy of solutions, while a poor one leads to nonconverging results, inaccurate solutions, and an increase in computation time. Due to the complexities of geometries in structural surfaces and the lack of automatic generation methods, it costs plenty of time to generate meshes in calculation processes. On the one hand, mesh quality is influenced by various factors, but on the other hand, it is unrealistic to adjust all of them to the best states in practice, which requires more computation power and computation time.

A good quality mesh should have a sufficiently large number of cells, implying that a further increase in the number of cells will not lead to significant changes in the solutions. The mesh needs to be fine enough in the regions where fluid characteristics change rapidly. Therefore, mesh properties need to be coordinated to make a trade-off between the accuracy and the time scale. In literature, different mesh properties, i.e., mesh types, domain sizes, and node distributions, are applied in RANS simulations of airfoil-shaped profiles. Lutton [5] indicated that the O-mesh is superior to the C-mesh in determining the pressure coefficients in the vicinity of leading and trailing edges, while the C-mesh shows a better resolution in the wake. Basha [6] constructed a structured O-mesh, a structured C- mesh, and a hybrid mesh on a NACA 0012 airfoil, while the hybrid mesh shows better results in drag coefficient predictions due to finer mesh resolutions in the boundary layer domain. Stuck et al. [7] selected a hybrid mesh with quadrilateral unstructured cells rather than triangular ones to achieve smoother transitions from the boundary mesh to the unstructured outer mesh for RANS simulations on a rudder profile. Parametric grid independence studies are conducted by varying chord-wise and layer-wise grids, respectively, and optimum values are chosen based on variations of integral force coefficients. Turnock et al. [8] carried out decoupled studies on the boundary location and mesh node distributions in different directions for a NACA 0012 section, while related parameters are determined by convergence study.

As no clear consensus about mesh properties was found in current studies, systematic parametric investigations on mesh properties are performed in this paper. As naval architects, our primary interest in CFD lies in applications related to the flow around ships and resultant hydrodynamic forces. In this particular case, we focus on the flow around the rudder and the resulting rudder performance. Kim et al. [9] stated that the NACA series are the most widely used economic rudder profiles. Thus, a classic NACA 0012 airfoil-shape rudder profile is utilized to study mesh properties for its simplicity in the shape, which makes mesh properties easier to change and enhances the applicability of obtained results for other complicated structures.

This paper analyzes the impacts of mesh properties on RANS simulations and proposes suitable meshing strategies for rudder hydrodynamics. Following this introduction, Section 2 further discusses the reason why a 2D rudder profile is selected as the research object in this paper. Section 3 explains definitions of mesh properties in this study. Section 4 illustrates applied research methods in detail, and Section 5 discusses the impacts of different mesh properties based on simulated cases. Finally, a recommended meshing strategy is presented to provide guidance for hydrodynamic analysis in Section 7.

#### **2. The Selection of the Research Object**

Bare hulls, propellers, rudders, or full-appendages ships were all studied in CFD methods in literature, though not all of them are suitable for mesh properties studies. There are great curvature changes on the surface of hulls and propellers, and that is why mixture mesh topologies and millions of cells are desired to capture their flow features. The complex nature makes it difficult to perform parametric investigations on a single mesh property economically. What's more, the geometries of hulls or propellers in the marine industry vary from each other significantly, and rules derived from limited cases may not be applicable for those with different configurations.

Although special rudders drew a growing interest, conventional spade rudders, or those modified but still geometrically similar versions like flap rudders, are still mainstream maneuvering devices, however. Besides, the meshing strategy aiming at a NACA profile can be easily applied to its geometrically similar counterparts. Force coefficients of rudders during ship maneuverings can be obtained based on 2D profiles in open water conditions, which also makes the investigation of 2D profiles beneficial to maneuverability studies. For all these reasons, in this paper, a 2D profile, in other words, a rudder with an infinite aspect ratio, is determined as the research object. According to Fujii [10], Fujii and Tsuda [11,12], the normal force coefficient of actual rudders whose aspect ratios are limited can be expressed as:

$$C\_N = 6.13 \sin \alpha\_R \frac{\Lambda\_G}{\Lambda\_G + 2.25} \tag{1}$$

where *CN* is the rudder normal force coefficient, and 6.13 is an approximate constant to estimate the partial derivative of *CN* to sin *<sup>α</sup><sup>R</sup>* of 2D sections. Liu et al. [2] calculated *<sup>∂</sup>CN ∂* sin *α<sup>R</sup>* for different rudder profiles by performing 2D CFD simulations, which can modify the original transformation relation in Equation (1) and obtain force coefficients for 3D rudders in open water conditions. Since rudders generally operate in the propeller slipstream, rudder hydrodynamics affected by the propeller can be estimated based on series model tests by Nienhuis [13], and the experimental data are shown in Figure 1. Compared with that of open water conditions, the stall angle for the rudder increased to about 35◦ with propeller impacts, while *<sup>∂</sup>CN ∂* sin *α<sup>R</sup>* was not significantly changed. Hence, based on 2D simulation results, corrected formulations for a rudder with an aspect ratio Λ*<sup>G</sup>* considering propeller impacts can be expressed as:

$$\begin{aligned} \mathsf{C}\_{L}^{\Lambda\_{G}} &= k\_{p} \left( \frac{\partial \mathsf{C}\_{L}^{2D}}{\partial \sin a\_{R}} \sin a\_{R} + \mathsf{C}\_{L\_{0}}^{2D} \right) \frac{\Lambda\_{G}}{\Lambda\_{G} + 2.25} \\\ \mathsf{C}\_{D}^{\Lambda\_{G}} &= k\_{p} \left( \frac{\partial \mathsf{C}\_{D}^{2D}}{\partial \sin a\_{R}} \sin a\_{R} + \mathsf{C}\_{D\_{0}}^{2D} \right) \frac{\Lambda\_{G}}{\Lambda\_{G} + 2.25} \end{aligned} \tag{2}$$

where *<sup>∂</sup>C*2*<sup>D</sup> L ∂* sin *α<sup>R</sup>* , *C*2*<sup>D</sup> <sup>L</sup>*<sup>0</sup> , *<sup>∂</sup>C*2*<sup>D</sup> D ∂* sin *α<sup>R</sup>* , *C*2*<sup>D</sup> <sup>D</sup>*<sup>0</sup> can be obtained from 2D CFD simulations, and *kp* can be estimated experimentally. The availability of Equation (2) in ship maneuvering studies was validated by Liu et al. [2,14] against two typical inland vessels. Apart from the transformation relation, another factor that makes the 2D study effective is that spanwise mesh generations can be omitted and cells in the profile section can be focused on. Therefore, 2D rudder profiles are set as the research object in following studies.

**Figure 1.** Comparison of the hydrodynamic coefficients of a rudder in open water conditions to those of the rudder in propeller slipstream. Data are adopted from Nienhuis [13].

#### **3. Mesh Properties**

In literature, mesh types and computation domain sizes are generally determined before studying grid independence on certain mesh parameters, and some settings are shown in Tables 1 and 2. Mesh properties in grid independence study are denoted as near-field and far-field parameters according to locations of distributions. All of these mesh properties are studied in this paper. The detailed descriptions of mesh properties above are in Sections 3.1–3.3.


**Table 1.** Mesh properties studies in literature: objects, Reynolds numbers, mesh types, and domain sizes.

**Table 2.** Mesh properties studies in literature: node distributions and turbulence models.


#### *3.1. Mesh Types*

Based on the connectivity of elements, meshes are identified as structured, unstructured, and hybrid. As regular connectivity can be expressed as a two/three-dimensional array, the structured model is highly efficient in the usage of computational resources. However, regular connectivity restricts the element type to quadrilateral in 2D and hexahedral in 3D, which inherently limits applications of the structured mesh for complex geometries. For simple geometries, a structured mesh may have better convergence and higher resolution than an unstructured mesh.

Instead of a single-block structured mesh, a block-structured or multiple-structured mesh is more widely applied. Three commonly used topologies, C-mesh, H-mesh, and Omesh, are shown in Figure 2. It is superior in computational efficiency to the unstructured mesh and more flexible in handling the complex geometries than the single-block structured

mesh. The choice of the mesh topology has a considerable influence on mesh quality [22], which depends on the domain geometry, the structure of the solution, and the topology in the adjacent block. According to Liseikin [22], in H-mesh, the computational domain is square and two singularities exist in the interior boundary. O-mesh and C-mesh are both in solid squares, while the interior boundary for the former one is smooth and the latter has a singularity the same type as the H-mesh. In this particular case, some modifications are made based on basic topologies above to accurately capture the characteristics of the flow near the rudder, as shown in Figure 2. For brevity, modified types are still marked as C-mesh, H-mesh, and O-mesh, respectively.

**Figure 2.** Original and modified mesh types for the block-structured mesh.

The irregular connection of unstructured meshes allows for more freedom in element choices, typically triangular in 2D and tetrahedral in 3D. Compared to that of the structured mesh, the unstructured mesh is more suitable for complex configurations, but it can be highly inefficient. Aftosmis et al. [23] reported that unstructured triangular meshes are 50 times more expensive in both memory and time than that of structured quadrilateral meshes with nearly the same accuracy. Unstructured meshes, however, have advantages over structured meshes in mesh refinement and generation time. In practice, hybrid meshes are applied instead of pure unstructured meshes. Hybrid meshes commonly consist of a structured portion near wall surfaces to capture the boundary layer and an unstructured portion that fills the domain. Thus, hybrid meshes inherit the advantages of both structured and unstructured meshes, such as good orthogonality to wall surfaces, suitability for mesh refinement, flexibility for complex geometries, and fast generation. Figure 3 demonstrates a hybrid triangle mesh, which has quadrilateral elements for the boundary layer and triangular elements for the remainder of the domain.

**Figure 3.** The mesh type of Hybrid-mesh.

#### *3.2. Domain Sizes*

The computational domain is the geometrical region that bounds the numerical simulation [17]. To obtain highly accurate solutions, the position of boundaries has to be discussed to demonstrate that the interior flow field is not disturbed. Thus, there is the need to study the domain size, which not only minimizes the influences of the boundaries, but also prevents an excessively large domain. Typical mesh types correspond to related topology structures, which lead to different shapes of calculation domains. The geometric parameters for different mesh types are shown in Figure 4, while subscripts, i.e., *in*, *out*, and *side*, stand for distances from the rudder to the inlet, outlet, and side boundary of corresponding mesh types, and *Or* is the radius of the domain of O-mesh.

#### *3.3. Node Distributions*

For structured meshes, since cell distributions vary with node distributions in different directions of the computation domain, convergence studies need to be done by varying related mesh parameters. In this study, node distributions are divided into two categories, including near- and far-field distributions.

Near-field distributions are located in the boundary layer regions near the wall, which are generally refined to capture the flow around the structure. The chord-wise spacing is the distance between two nodes along with the profile, and the layer-wise growth rate is the distances between adjacent mesh points along a mesh line, which determine the boundary layer mesh along with the first mesh height.

Different mesh types share similar near-wall distributions, while it is difficult to define far-field distributions owing to differences in mesh topologies. In this paper, far-field distributions are defined as node distributions that do not contain near-field distributions.

**Figure 4.** The geometric parameters for different mesh types.

#### **4. Research Approach**

A series of parametric studies are conducted to observe the impacts of mesh properties on the prediction of the accuracy of solutions. The classic validation profile NACA 0012 is tested through commercial CFD code ANSYS Fluent. Strategies for investigating mesh properties mentioned in Section 4.1. Numerical methods utilized to study the hydrodynamic performance of the rudder profile are described in Section 4.2, and the error calculation method is mentioned in Section 4.3.

#### *4.1. Strategies*

Different mesh properties contribute to mesh quality in varying degrees, and it is the combination of these properties that determines the accuracy of calculation results. Instead of rough global refinement of cells [24], several main impact factors on mesh fineness and the required number of cells are illustrated with simulations. This parametric mesh independence study is chosen because optimization on individual parameters specifically addresses the crucial impact factor of mesh quality, avoiding the waste of cells [7]. The global refinement is more frequently applied for unstructured meshes than structured meshes. However, due to the lack of specific classification on mesh properties, mesh properties above are not independent of each other strictly. Cases with different mesh types also differ from each other in domain sizes, node distributions, and numbers of mesh elements, etc.

As shown in Figure 5, the typical mesh type corresponds to a related topology structure, which leads to different shapes of calculation domains and connection patterns between the profile and domain boundaries. With domain size varying, node distribution in different directions changes. Further, the number of mesh elements is concerned with all mesh properties above.

**Figure 5.** Relationship between different mesh properties.

The strategy in Figure 6 is applied to investigate how different mesh properties affect the calculation results of rudder hydrodynamics. Firstly, cases of C-mesh type with Reynolds numbers in a wide range are tested to obtain the critical Reynolds number (*Recritical*), above which hydrodynamic coefficients tend to be stable. Secondly, the domain sizes of the C-mesh are varied to obtain the proper boundary sizes that are large enough to eliminate boundary effects for the following simulations. Thirdly, cases are divided into 4 groups classified by 4 different mesh types. Since it is inappropriate to compare mesh qualities of different mesh types when other mesh properties are not determined yet, the impacts of node distribution on each mesh type should be studied independently. In this part, the node distribution is the only factor that affects mesh quality for the typical mesh type, and the effects of node distributions are studied based on *Recritical* and the domain size obtained in the previous stage. Therefore, grid independence studies are conducted on parameters in the near-field, i.e., chord-wise spacings and layer-wise growth rates, and the far-field, i.e., radial nodes, wake nodes, and element sizes. Finally, four mesh types with optimal mesh properties are compared to give recommended mesh settings for rudder hydrodynamics.

#### *4.2. Numerical Methods*

The RANS methods presented in this paper use a time-average Reynolds decomposition, which assumes that all the components of the flow velocity and pressure consist of a mean value and a bounded fluctuation to represent turbulence. The governing equations for incompressible viscous flow are as follows:

$$\frac{\partial \bar{u}\_i}{\partial x\_i} = 0 \tag{3}$$

$$
\rho \frac{\partial \overline{u}\_i}{\partial t} + \rho \overline{u}\_j \frac{\partial \overline{u}\_i}{\partial x\_j} = f\_i - \frac{\partial \overline{p}}{\partial x\_i} + \frac{\partial}{\partial x\_j} \left( \mu \frac{\partial \overline{u}\_i}{\partial x\_j} - \rho \overline{u\_i' u\_j'} \right), \quad (i = 1, 2, 3) \tag{4}
$$

where *u*¯*<sup>i</sup>* is the mean velocity component in the *x* direction, *ρ* and *μ* are the density and viscosity coefficient of the fluid, *P*¯ is the mean pressure, *fi* is the body force, while *ρu i u j* is the Reynolds stress. In RANS simulations, turbulence models are introduced to close equations above. The *k*-*ω* SST model was designed to give highly accurate predictions of the onset and the amount of flow separation, which provides a better description of flow around the rudder. Therefore, the *k*-*ω* SST model is applied for simulations in the following parts. Numerical settings for ANSYS Fluent is shown in Table 3.

**Figure 6.** Strategy applied to study different mesh properties.

**Table 3.** Numerical settings for RANS simulations with ANSYS.


#### *4.3. Error Calculation*

Routinely, nondimensional coefficients are used to compare the rudder hydrodynamic performance with that of various design choices, while main hydrodynamic characteristics are the lift coefficient *CL*, the drag coefficient *CD*, and the moment coefficient *CM* ([1]), which are expressed as:

$$\mathbf{C}\_{\rm L} = L\_{\rm R} / \left( 0.5 \rho V\_{\rm R}^2 A\_{\rm R} \right) \tag{5}$$

$$\mathbf{C}\_{\rm D} = D\_{\rm R} / \left( 0.5 \rho V\_{\rm R}^2 A\_{\rm R} \right) \tag{6}$$

$$\mathcal{L}\_{\mathrm{M}} = M\_{\mathrm{R}} / \left( 0.5 \rho V\_{\mathrm{R}}^2 A\_{\mathrm{R}} c \right) \tag{7}$$

where *L*R, *D*<sup>R</sup> and *M*<sup>R</sup> are the lift force, drag force, and moment, respectively, *ρ* is water density, *α* is the angle of attack, *V*<sup>R</sup> is the rudder inflow speed, *A*<sup>R</sup> is the rudder area, and *c* is the chord length of the rudder profile. The center of pressure is set as 0.25 *c* from the rudder's leading edge, based on engineering practice and experimental data [25].

To evaluate the mesh quality, wind tunnel results of the NACA 0012 profile, which are tested by Ladson [25], are taken as the benchmark case. Moreover, three independent compressible CFD codes [19], i.e., CFL3D (NASA LaRC, USA), FUN3D (NASA LaRC, USA), and NTS (NTS, Russia), are cited to make a crosswise comparison with the Fluent solutions. Results of the benchmark and the CFD codes are achieved in essentially incompressible air with a Mach number of 0.15. Thus, they are comparable with the results obtained in incompressible water in this paper. Presented results are compared to that of the benchmark cases in relative differences under various angles of attack at a critical Reynolds number as follows:

$$
\Delta \mathbb{C}\_a = \frac{\mathbb{C}\_{aFunc}}{\mathbb{C}\_{aBernular}} - 1 \tag{8}
$$

$$\Delta \text{C(AVE)} = \frac{\sum\_{i=1}^{n} |\Delta \mathcal{C}\_{a\_i}|}{n}, \quad i = 1, 2, \dots, n \tag{9}$$

where Δ*C<sup>α</sup>* is the relative differences of the calculated lift, drag, or moment coefficients achieved by Fluent and provided by the benchmark data, respectively, and Δ*C*(AVE) is the average relative differences for a series of *α*.

#### *4.4. Uncertainty Estimation*

The uncertainty estimation is important to evaluate the reliability of data in CFD calculations. According to International Towing Tank Conference [26], simulation numerical uncertainty *USN* and simulation modeling uncertainty *USM* need to be assessed respectively. In this paper, only the numerical errors are considered since modeling errors are difficult to quantify in practical applications [27]. Generally, numerical errors relate to grid size, time step, and other parameters. This paper performs steady simulations and focuses on mesh properties. Accordingly, only the grid uncertainty *UG* is calculated for uncertainty estimations. With validation uncertainty *UV*, which is simplified as *UG* in this study, obtained, *UV* and the comparison error |*E*| can be compared. If |*E*| < *UV*, the validation is achieved at *UV* level.

The Grid Convergence Index (GCI) method [28,29] and the correction factor (CF) based method [30] remain popular in uncertainty estimations, however, it is deficient that both methods fail to give reliable estimations when the order of grid convergence *p* is too high (i.e., *p* > 2) or too small (i.e., *p* < 0.5). In other words, GCI and CF methods require data in the so-called asymptotic range, which means that data on grids are fine enough to give a single dominant term in a power series expansion of the error [31]. Considering that practical flow phenomenons are generally complex, and it is inevitable that *p* may be not in the range of [0.5, 2], Eça and Hoekstra [31] proposed another GCI method based on least-square roots (LSR-GCI). Such a method allows for error estimations even if the monotonic convergence is not achieved, and it is proven to be more reliable compared with GCI and CF based methods by De Luca et al. [32]. Therefore, the LSR-GCI method [31] is applied for uncertainty estimations in this paper.

For four solutions *φi*,(*i* = 1, 2, 3, 4) obtained from systematically refined cases and characterized by typical cell sizes *hi*,(*i* = 1, 2, 3, 4), the discretization error can be estimated with one of the four following equations:

$$
\varepsilon\_{\phi} \simeq \phi\_i - \phi\_0 = \alpha h\_i^p \tag{10}
$$

$$
\varepsilon\_{\Phi} \simeq \phi\_{\bar{i}} - \phi\_0 = ah\_{\bar{i}} \tag{11}
$$

$$
\varepsilon\_{\phi} \simeq \phi\_i - \phi\_0 = \alpha h\_i^2 \tag{12}
$$

$$
\varepsilon\_{\Phi} \simeq \phi\_{\bar{i}} - \phi\_0 = \kappa\_1 h\_{\bar{i}} + \kappa\_2 h\_{\bar{i}}^2 \tag{13}
$$

where *φ*<sup>0</sup> is the estimate of the exact solution, *α* is an undetermined constant, and *p* is the observed order of grid convergence. Besides, the following functions are proposed to select the error estimator:

$$S\_{RE}(\phi\_0, \alpha, p) = \sqrt{\sum\_{i=1}^{n\_{\mathcal{S}}} \left(\phi\_i - \left(\phi\_0 + ah\_i^p\right)\right)^2} \tag{14}$$

$$S\_{\rm RE}^w(\phi\_0, a, p) = \sqrt{\sum\_{i=1}^{n\_{\rm f}} w\_i \left(\phi\_i - \left(\phi\_0 + a h\_i^p\right)\right)^2}, (w\_i = \frac{\frac{1}{h\_i}}{\sum\_{i=1}^{n\_{\rm f}} \frac{1}{h\_i}}) \tag{15}$$

$$S\_1(\phi\_0, \boldsymbol{\alpha}) = \sqrt{\sum\_{i=1}^{n\_{\mathcal{S}}} (\phi\_i - (\phi\_0 + ah\_i))^2} \tag{16}$$

$$S\_1^w(\phi\_0, a) = \sqrt{\sum\_{i=1}^{n\_\mathcal{S}} w\_i (\phi\_i - (\phi\_0 + a h\_i))^2}, (w\_i = \frac{\frac{1}{h\_i}}{\sum\_{i=1}^{n\_\mathcal{S}} \frac{1}{h\_i}}) \tag{17}$$

$$S\_2(\phi\_0, a) = \sqrt{\sum\_{i=1}^{n\_{\mathcal{S}}} \left(\phi\_i - \left(\phi\_0 + ah\_i^2\right)\right)^2} \tag{18}$$

$$\mathcal{S}\_2^w(\phi\_0, a) = \sqrt{\sum\_{i=1}^{n\_\mathcal{S}} w\_i \left(\phi\_i - \left(\phi\_0 + a h\_i^2\right)\right)^2}, (w\_i = \frac{\frac{1}{h\_i}}{\sum\_{i=1}^{n\_\mathcal{S}} \frac{1}{h\_i}}) \tag{19}$$

$$S\_{12}(\phi\_0, \alpha\_1, \alpha\_2) = \sqrt{\sum\_{i=1}^{n\_\mathcal{S}} \left(\phi\_i - \left(\phi\_0 + \alpha\_1 h\_i + \alpha\_2 h\_i^2\right)\right)^2} \tag{20}$$

$$S\_{12}^{w}(\phi\_0, a\_1, a\_2) = \sqrt{\sum\_{i=1}^{n\_{\mathcal{S}}} w\_i \left(\phi\_i - \left(\phi\_0 + a\_1 h\_i + a\_2 h\_i^2\right)\right)^2}, (w\_i = \frac{\frac{1}{h\_i}}{\sum\_{i=1}^{n\_{\mathcal{S}}} \frac{1}{h\_i}}) \tag{21}$$

Details in selecting the error estimator are shown in Figure 7 according to Eça and Hoekstra [31]. *φ*∗ 0,*RE*, *α*<sup>∗</sup> *RE*, *p*<sup>∗</sup> *RE* are variables that make *SRE*(*φ*0, *α*, *p*) minimum, which can be denoted as  *φ*∗ 0,*RE*, *α*<sup>∗</sup> *RE*, *p*<sup>∗</sup> *RE* = *S*−<sup>1</sup> *RE* min{*SRE*(*φ*0, *α*, *p*)} with the inverse function symbol *S*−<sup>1</sup> *RE*. Similarly, we have *<sup>φ</sup>w*<sup>∗</sup> 0,*RE*, *<sup>α</sup>w*<sup>∗</sup> *RE*, *<sup>p</sup>w*<sup>∗</sup> *RE*, *φ*<sup>∗</sup> 0,1, *α*<sup>∗</sup> <sup>1</sup>, *p*<sup>∗</sup> <sup>1</sup>, *<sup>φ</sup>w*<sup>∗</sup> 0,1, *<sup>α</sup>w*<sup>∗</sup> <sup>1</sup> , *<sup>p</sup>w*<sup>∗</sup> <sup>1</sup> , *φ*<sup>∗</sup> 0,2, *α*<sup>∗</sup> <sup>2</sup>, *p*<sup>∗</sup> <sup>2</sup>, *<sup>φ</sup>w*<sup>∗</sup> 0,2, *<sup>α</sup>w*<sup>∗</sup> <sup>2</sup> , *<sup>p</sup>w*<sup>∗</sup> 2 , *φ*∗ 0,12, *α*<sup>∗</sup> <sup>12</sup>, *p*<sup>∗</sup> <sup>12</sup>, *<sup>φ</sup>w*<sup>∗</sup> 0,12, *<sup>α</sup>w*<sup>∗</sup> <sup>12</sup> , *<sup>p</sup>w*<sup>∗</sup> <sup>12</sup> . Corresponding standard deviations can be expressed as:

$$\sigma\_{RE} = \sqrt{\frac{\sum\_{i=1}^{n\_{\mathcal{S}}} n w\_i \left(\phi\_i - \left(\phi\_0 + a h\_i^{p'}\right)\right)^2}{\left(n\_{\mathcal{S}} - 3\right)}}\tag{22}$$

$$\sigma\_1 = \sqrt{\frac{\frac{\sum\_{i=1}^{n\_{\mathcal{S}}} n w\_i \left(\phi\_i - \left(\phi\_0 + a h\_i\right)\right)^2}{\left(n\_{\mathcal{S}} - 2\right)}}\tag{23}$$

$$\sigma\_2 = \sqrt{\frac{\frac{n\_{\mathcal{S}}}{\sum\_{i=1}^{n\_{\mathcal{S}}} m v\_i \left(\phi\_i - \left(\phi\_0 + a h\_i^2\right)\right)^2}{\left(n\_{\mathcal{S}} - 2\right)}}\tag{24}$$

$$\sigma\_{12} = \sqrt{\frac{\sum\_{i=1}^{n\_{\mathcal{S}}} n w\_i \left(\phi\_i - \left(\phi\_0 + \alpha\_1 h\_i + \alpha\_2 h\_i^2\right)\right)^2}{\left(n\_{\mathcal{S}} - 3\right)}}\tag{25}$$

For no-weighted cases, *nwi* = 1. Following criteria are applied in selecting the error estimators:


**Figure 7.** Procedure for error and uncertainty estimations by GCI-LSR method recommended by Eça and Hoekstra [31], De Luca et al. [32].
