Contaminant removal efficiency of floating treatment wetlands

Floating treatment wetlands are new ecological infrastructures for stormwater treatment. Despite a recent proliferation in their usage, their contaminant removal efficiency e continues to draw research attention. Here, the e from idealized FTWs is numerically computed across a wide range of flow and geometric conditions while accommodating joint contributions of advection, turbulent dispersion, and vegetation removal. The emerging mathematical structure describing e bears resemblance to a simplified plug flow model and supports an empirical shallow-basin model from long-term field measurements. The present model indicates that e remains significantly influenced by a Dämkohler number that quantifies the effects of both vegetation and flow properties. The impacts on e of the underflow region and contaminant blockage on the removal mechanisms are also investigated.


Introduction
Floating treatment wetlands or islands (FTWs or FTIs) are recent additions to a proliferating ecological infrastructure for stormwater and wastewater treatment (Van de Moortel et al 2010, de Stefani et al 2011, Braskerud 2002, Garcia et al 2010, Haberl et al 2003, Langergraber and Simunek 2005, Holland et al 2004, Dierberg et al 2005, Keefe et al 2004, Kumar and Zhao 2011, Afzal et al 2019, Colares et al 2020. FTWs, schematically featured in figure 1, consist of an underflow layer and a porous canopy layer where much of the actual contaminant removal occurs. The underflow layer is usually turbulent and impacted by vegetated elements (mainly roots). The underflow layer is occasionally referred to as an inner layer because it is bounded by two shear-stresses: a wall stress (at the ground) and a stress at the vegetation-no vegetation interface (at the upper end). The premise of many FTWs rests on the fact that the canopy layer may be approximated by a porous medium (Plew 2011, Huai et al 2012, Li et al 2019 where the flow is sufficiently slow (but can be turbulent) thus allowing for removal mechanisms and chemical transformations to occur. However, the presence of a turbulent underflow layer complicates the evaluation of contaminant removal efficiency (e) by FTWs, the compass of this letter. Velocities in the underflow region are usually high and can 'short-circuit' the vegetation removal mechanism of contaminants. The short-circuiting refers to contaminants traversing a path in the FTW without encountering a removal mechanism (Kadlec 2000), which is assumed to be the vegetation sink term here. On the other hand, this short-circuiting mechanism is ameliorated by turbulent dispersion within the underflow region. An enhanced eddy diffusivity (D t ) allows for contaminants to be vertically (and efficiently) transported from the underflow region into the overlying vegetation layer where vegetation removal mechanisms can operate. This turbulent mass transport occurs along a mean vertical concentration gradient. Moreover, the turbulent state of the underflow region leads to higher D t and shear stresses within the vegetated layer thereby enhancing the vegetation removal efficiency itself. The scope of this work is to explore these competing effects arising from the presence of an underflow region on e in FTWs. Even within this restricted scope, many assumptions and simplifications must be made about the transport processes and the removal mechanisms impacting e. For the latter, two end-member cases where the removal mechanism experiences no blocking and appreciable blocking due to the gradual accumulation of contaminants around the removal sites are considered. Also, simplified D t models and their associated mixing length characterization are employed to link the underflow region to the vegetated zone, where contaminant removal ultimately occurs.
When the dominant terms in the contaminant mass balance are advective transport (plug flow or infinite Peclet number) and a first-order removal mechanism, the balance reduces to where k v is an inverse removal time scale, U p is a mean plug flow velocity, C is the time-averaged contaminant concentration, and x is the longitudinal direction. Throughout, all the variables use the International System of Units (SI) unless otherwise stated. The C variations are given by which leads to a removal efficiency (Carleton 2002, Kadlec and Wallace 2008, Khan et al 2013) where C(0) and C(L) are the concentrations at the FTW entrance (x = 0) and outlet (x = L), and Da = k v (L/U p ) is a D amkohler number, L is the FTW length, and L/U p is considered the appropriate advective time scale. For large L and slow U p , Da becomes large and the removal efficiency is enhanced for a preset k v . This representation of e is consistent (at least in mathematical form) with a large corpus of experiments and detailed models (Kadlec and Wallace 2008). Other representations emphasize the role of aspect ratio (L to the width W) instead of Da and propose (Thackston et al 1987) This expression was derived from data fitting to longterm field experiments collected in a variety of large shallow wetlands of diverse sizes, shapes and vegetation type. The expression suggests that for long and narrow wetlands (i.e. large L/W), the maximum operationally achievable e is 0.85. It also highlights the complication of using equations (1) to (3) to infer e when the flow field is highly simplified. In steadyuniform flow where the volumetric inflow rate Q in per unit W is held constant and Manning's formula is used to estimate U p for large W yields U p ∝ W −2/5 (see appendix). Equating the exponential terms in equations (3) and (4) require U p to scale as W +1 . Such analysis cannot even predict the scaling exponent between U p and W suggested by the experiments summarized in equation (4). A tenuous connection between Manning's roughness n f and geometry (i.e. n f scales as W -7/4 necessitating a water depth scaling as W 9/20 ) may recover U p scaling as W +1 but this connection has no clear rationale (see appendix). For FTWs, contaminant removal efficiency encounters another level of complexity partly because of the presence of a turbulent underflow region. The goal is to develop a parsimonious mathematical representation for e that is commensurate in complexity to equation (3) but informed by detailed model calculations that accommodate many (though not all) processes that are common to all FTWs (including the underflow zone). This work is by no means sufficient to address all issues associated with FTW operation or design; however, it is a necessary first step to advance understanding of contaminant removal efficiency in FTWs.
As a starting point, the simplified FTW sketched in figure 1 is considered. Much of the description of the flow for this idealized FTW has been explored elsewhere (Li et al 2019) and are not repeated in detail. Briefly, the FTW system consists of two layers: the canopy layer and the underflow region. Contaminants enter at x = 0 uniformly in z. They are transported and removed due to the mutual effects of advection, turbulent dispersion (ignored in equation (1)), and vegetation removal (other biogeochemical transformations or removal at the ground can be accommodated here but are momentarily not considered). Based on the aforementioned two-layer configuration, a model for the description of contaminant removal efficiency is now sought. Numerical model runs on contaminant removal processes are conducted by varying the FTW dimensions and the vegetation removal mechanisms through the sink term. The basic premise of these runs is to ensure sufficient variability in the key dimensionless numbers relevant to contaminant transport removal, namely the Peclet number, the D amkohler number, and geometric properties such as relative length scales describing the FTW configuration.

Theory
To estimate the contaminant removal efficiency, the mass transport equation that includes advective and dispersive transport as well as removal (S r ) of contaminants by the vegetation layer is considered. When S r represents chemical transformations, then it can be finite in all layers. Here, S r is purposely restricted to the vegetation layer where the uptake may be due to physical (e.g. deposition), chemical (e.g. sorption onto the root surface), or biological (e.g. via passive or active plant uptake or via removal in the biofiltration medium and roots) depending on the contaminants. These contaminants commonly include total nitrogen (TN), total phosphorus (TP), total suspended solids (TSS), metals and so on. For simplicity, steady-state conditions are assumed and the inflow contaminant C(0) is set to be vertically uniform so as to facilitate the determination of e across different numerical runs that vary in FTW sizes and key dimensionless numbers. Also, such assumptions facilitate comparisons with equation (3). For a large wetland width, symmetry considerations in the lateral direction permit a 2D simplification so that the conservation of contaminant mass becomes, where t is time, F x and F z are the longitudinal and vertical mass fluxes respectively, u(z) is the local mean velocity in the streamwise direction, D is an apparent dispersion coefficient representing the summed effects of molecular (D m ) and D t that requires a mathematical model ( The treatment of scalar transport as 2D whereas the treatment of u(z) as 1D (via the measurements) may raise objections. However, a number of studies have already shown that the momentum adjustment distance (i.e. the distance where ∂u/∂x ≈ 0 after x = 0) is much shorter (by factors of 5-10) than its scalar counterpart (Siqueira et al 2012). For a 1D incompressible flow (i.e. ∂u/∂x = 0), the continuity equation (∂u/∂x + ∂w/∂z = 0) leads to F z = −D∂C/∂z since the vertical velocity (w) must be zero at the FTW floor (w = 0 at . The vegetation removal term, S r , is assumed to be a piece-wise function, which is set to zero for the underflow region. This assumption ignores other physical (e.g. deposition and removal at the FTW floor) and chemical transformations within the underflow region relative to the canopy layer. In the canopy layer, S r is modeled using two different forms, where C eq (>C o ) is a contaminant concentration when the vegetation removal is suppressed (i.e. at C eq , S r is maximal for r = 1; S r = 0 for r = 2), k v is an inverse time scale that encodes all the physical and biochemical processes occurring in the vegetation layer (Chang et al 2012, Chua et al 2012, Xavier et al 2018, and C o is the concentration where S r is maximum for r = 2. When r = 1, the vegetation removal is a first-order process common to many contaminant removal studies at low concentrations ( . However, as the C becomes high, first-order models are expected to overestimate the removal rate. Here, we term this phenomena as blockage and it means that the removal rate decreases as C becomes large because of competition between the many contaminant molecules for the removal surface area. That is, the presence of nearby molecules near the sink 'blocks' access of other molecules further away. A second-order model (parabolic, when r = 2 in equation (6) so as to resemble competition in a logistic equation) is employed to accommodate such an effect as shown in figure 1. The (1 − C/C eq ) r−1 (referred to as ϕ c with ϕ c < 1) in the r = 2 removal term can be understood as a reduction factor to k v , hence implying an intrinsic inverse time scale of (k v ϕ c ). The goal here is to contrast two differing formulations for S r (r = 1, 2) and track this contrast all the way to e. The case where r = 0 (i.e. a constant S r independent of C) is not considered but can be accommodated in the formulation here. That is, in analogy to reaction kinetics, r measures the order of the reaction rate describing S r . It is noted that wetland removal is a complex process that includes many factors e.g. microbially mediated processes, chemical and nutrient removal processes, sedimentation, solute sorption, passive and active plant uptake among many others (Polprasert and Khatiwada 1998, Reddy and D'angelo 1997, Kadlec and Wallace 2008, Singh et al 2019. To assess generic formulation for the removal efficiency, it is convenient to lump these processes into a single expression for S r encoded by k v . Needless to say, there are many instances where increased complexity in the removal process prevents a single k v to be defined. There are also instances where increasing the complexity of the S r model need not yield improved predictive skills (Paudel and Jawitz 2012). The S r formulations may need revisions for a problem-specific case study. In the case of biofilm removal modeling, k v can be interpreted as a first-order reaction rate describing how dissolved materials in the vicinity of solid surface are transported into the biofilm (Kadlec and Wallace 2008). Clearly, in this case, k v pertains to the physical properties of water and biofilm (e.g. diffusion coefficient and biofilm thickness) and the removal sites on stems and roots (e.g. canopy density, and surface area) where the biofilm forms and grows (Khatiwada andPolprasert 1999, Kadlec andWallace 2008). Chemically, the removal process in the wetlands can involve more than one reaction and more than one species so that k v must be understood as an effective first-order reaction rate parameter. The S r formulation can also be applied to model sorption mechanism using partition coefficients (k p 's), e.g. in a linear sorption model, C s = k p C; in Freundlich model, C s = k p1 C p ; and in Langmuir model C s = k p2 ( C C+b ), where C s is the sorption term (can be equivalent to S r ), p and b are two coefficients representing the sorptive process (Kadlec and Wallace 2008, Clark 2011). Finally, k v can also be interpreted as the nutrient uptake inverse time scale (Kadlec and Wallace 2008) when modeling nutrient uptake by plants since the dissolved material exchange is most intense near the root region.
For the apparent dispersion coefficient, it is assumed that D = D t + D m , where the molecular diffusion in water D m is much smaller than the fluid viscosity (molecular Schmidt number, Sc, far exceeds unity). However, the turbulent Sc is of order unity such that the turbulent diffusivity D t = ν t /Sc ≈ ν t . Many studies have reported some deviations from unity in turbulent Sc, especially in the roughness sublayer of vegetation (Flesch et where the selection of l m (z) in the underflow and the vegetation zones is discussed in the appendix. Finally, for the purposes of selecting model runs so as to vary the relative importance of scalar transport and removal mechanisms, equation (5) is expressed in dimensionless form, where C ′ = C/C o is the normalized concentration,  (8) suggests that the removal efficiency of the FTWs be expressed as a generic function Here, it is assumed that D ′ scales with bulk flow variables similar to those already included in equation (9). An expression for f (.) similar to equation (3) is now sought at steady-state reflecting the processes in equation (8). By construct, f (.) cannot include all the complex processes (biotic and abiotic impacting S r ) that might be active in FTW outside those considered in equation (8). However, it is safe to state that the processes accommodated in equation (8) are, at minimum, required in evaluating the efficiency of all FTWs.

Results
The FTW geometry and flow velocity for the numerical runs were chosen to match the u(z) measured in laboratory experiments (Plew 2011). This choice ensures that the mean flow field and turbulent diffusivities driving contaminant transport are plausible and commensurate with detailed measurements. In these experiments, the flow depth is set as H = 0.2m and the relative thickness of the flow is described by For the FTW geometry, 13 different wetland lengths (different Pe L ) respectively denoted by ε −1 = 1, 3, 5, 7 9, 10, 12, 14, 16, 18, 20, 22 and 25 are investigated. In all these cases, it is assumed that the measured u(z) is not impacted by FTW edges and maintains its 1-D character (for simplicity). The C eq is set arbitrarily to C eq = 2 C| x=0 since the canopy type and density are fixed in the laboratory experiments. In total, 832 runs were conducted and the outcome of these runs are used to analyze a reducedorder model for e (i.e. determining f (.)). The numerical procedure and results for e are featured in the appendix. For evaluation of e, steady-state conditions (∂C ′ /∂t ′ = 0) are presumed as before when solving for C. Two representative concentration contours of C at steady-state for r = 1 and r = 2 are presented in figure 2 to illustrate the role of blocking. Figure 2 shows how C decreases from inlet to outlet after attainment of steady-state. Only the canopy zone is allowed to remove contaminants in the model calculations. The underflow region plays a dual role: allowing contaminants to advect longitudinally (i.e. short-circuiting) but higher l m and turbulent diffusivities (see appendix) allow contaminants to be vertically transported into the canopy region. Once in the canopy zone, the removal mechanism acts to remove contaminants there. Contaminant removal increases with increasing k v as expected from equation (3). Moreover, when different removal kinetics (r) are explored, the linear (r = 1) mechanisms are more efficient in removing contaminants when compared to the case with blocking (r = 2) (discussed later). Finally, sample contaminant concentration contours with different FTW length (L) are shown in figure 3 for illustration. Table 1. Summary of the reported krp for wetlands (not necessarily FTW) in the literature. It is noted that the literature inverse time scales reported are based on the time interval between two consecutive concentration measurements assuming well-mixed conditions in the entire wetland volume V. This implies that krp describes a system removal rate that can be related to kv using kv = krpV/Vc, where V is the entire volume of the wetlands and Vc is the volume of the active layer removing mass (in the FTW does not include the underflow region). Therefore, we set 10 −5 as the smallest value for kv in the numerical solutions to reflect this volume difference. It is found later in the numerical runs that values smaller than 10 −5 s −1 have negligible impact on the computed e.

Author
System type Material krp (s −1 ) e (%) Chua et al ( Figure 3 demonstrates that the contaminant concentration in the canopy layer is much lower than that in the underflow layer setting a vertical gradient for mass transport directed from the underflow into the canopy region. By comparing different wetland length (L), it is evident that when L increases, the removal efficiency also increases (partly by allowing for vertical transport to operate). However, it is less significant when compared with increases in k v , as shown in figure 2.
After examining correlations between e and different parameters in equation (9), a concise model for contaminant removal efficiency emerges as (details in appendix) where c 1 and c 2 are two empirical coefficients related to FTW configurations and canopy properties (e.g. ε g , a). The coefficient c 1 (≤1) can be interpreted as a maximum background efficiency that can be attained for such flow-vegetation-geometry configuration, whereas c 2 can be interpreted as the necessary amplification to Da L needed to match the plug-flow formulation when c 1 = 1 (so as to compensate for the short-circuiting). The Da L in equation (10) accommodates the joint contribution of the vegetation removal and longitudinal transport. For clarity, a linearly equivalent removal kinetic coefficient k * v instead of k v is utilized in formulating Da L for r = 2. The parameter k * v = k v (1/3 + 2/3ϕ c ) can be understood as an equivalent linear slope of the described second-order model derived from the integral constraint´C 0 k * v CdC =´C 0 S 2 dC and the terminal condition C = C eq . The two coefficients c 1 and c 2 , describing the deviation of the FTW efficiency from equation (1) arise due to the interplay between the vegetation removal mechanism and the underflow layer. As earlier noted, c 1 and c 2 must reflect the FTWs configuration and short-circuiting impact due to contaminants entangled in the underflow layer bypassing the vegetation removal mechanism. This implies that these coefficients vary with flow thickness ε g and canopy density a. The comparison between the fully modeled and the fitted f (.) is shown in figure 4 for all the runs. Figure 4 demonstrates that FTW efficiency can be reasonably described by the simplified expression featured in equation (10), implying the removal process is dominated by the D amkohler number (i.e.  (10) are shown in black solid and dashed lines for r = 1 and r = 2, respectively. It is noted here for r = 2, DaL is calculated based on k * v . Right: The overall fit of the predictive model. The inset is a zoomed-in to emphasize agreement for small e. The data are labeled with the same markers as the left and the colorbar indicates data density. The one-to-one line is also shown for reference.
Da L = k v U/L for the r = 1 case, or Da L = k * v U/L for the r = 2 case). The mathematical form of the simplified expression agrees with the plug flow case in equation (3) albeit with different coefficients. For modeling purposes, c 1 and c 2 here are temporarily optimised by fitting the predictive model to the numerical model runs. Analysis (and fitted values in the appendix) show that c 1 and c 2 are related to ε g (as expected for FTW). The r = 1 model is about 30% more efficient in terms of removing contaminants compared with that of r = 2 (with blockage) when the D amkohler number is small (i.e. Da L = k v U/L < 0.5). For large D amkohler numbers, the difference between the two vegetation removal models represented by r = 1 and r = 2 are significant. As before, high removal kinetics (k v or k * v ) and long FTW length boost contaminant removal as anticipated by equation (3). It is worth noting that for small ε g (and r = 1), c 1 ≈ 1 and c 2 ≈ 1 consistent with plug-flow model calculations (see appendix). Also, the maximum empirical efficiency of 0.85 approximated from figure 4 is reasonably bounded by c 1 (and r = 1) when ϵ g < 0.25. Maximal e values computed here for lower k v and r = 1 (e varies from 0.55-0.75) are commensurate with the higher-end measured values reported in table 1 lending some confidence to the plausibility of the model.

Conclusion
Contaminant removal efficiency in FTWs is analyzed and a concise model that collapses numerical results is proposed. The model indicates that the contaminant removal efficiency in FTWs remains primarily controlled by D amkohler number (as with conventional wetlands) but with some modifications. These modifications account for the joint influence of FTWs configuration and vegetation removal effects through two coefficients c 1 and c 2 , which accommodate the interplay between the removal mechanism and the underflow layer. Different from the plug flow case where a 100% removal is theoretically attainable, the present model indicates that the removal efficiency of FTWs cannot exceed the background efficiency (c 1 ), supporting prior findings (Thackston et al 1987, Afzal et al 2019. The model also shows that the removal efficiency is boosted as k v increases, which is consistent with the results from TP measurements summarized in table 1. Furthermore, the numerical runs here demonstrate that the contaminant blockage can decrease the removal efficiency by ameliorating the vegetation removal mechanisms. While the model results remain confined to a naive representation of FTWs, they highlight the roles of the underflow region (and factors that ameliorate its short circuiting effect) and blockage in reducing contaminant removal efficiency. Finally, the model runs can be used to assist in the designs of experiments targeting the removal efficiency of FTWs or hypotheses about new mechanisms governing their efficiency.

Data Availability
All data used for the flow conditions were published in the literature. All model results of the scalar concentration analysis are available at https://doi.org/10.7924/r48k7bv5t.

Acknowledgments
S L thanks R Callari-Kaczmarczyk for providing editorial comments on a prior version of this work and acknowledges the financial support from the Nicholas School of the Environment at Duke University. G K acknowledges support from the U.S. National Science Foundation (NSF-AGS-1644382, NSF-IOS-1754893, and NSF-AGS-2028633 and When S o and n f are constants, such analysis yields U p ∝ W −2/5 and y d = W −3/5 . The analysis here also shows that a U p ∝ W +1 can only be recovered when n f ∝ W −7/4 resulting in y d ∝ W 9/20 .  respectively. We conducted a separate analysis to determine the optimal values of the three proportionality coefficients. In this analysis, the modeled  shear stress profiles derived from measured u(z) and the linear-mixing length model are compared to those reported by Plew (2011). To further reduce the degrees of freedom in parameter determination, the k IV and k III are assumed to be the parameters to be optimized but constrained by the ranges provided by Huai et al (2012). The parameter k II is separately determined by continuity considerations in l m within zones III and zone I (and thus no a priori constraints are enforced). The optimized values for these coefficients are reported in table B1 for the differing geometry. As shown from table B1, k II (determined from l m continuity considerations only) is generally near the expected range [0.2, 0.4], which is not a priori enforced. The value of k IV is close to the expected von Karman constant ( ≈ 0.4) for all ε g runs thereby lending some confidence in the optimized parameter selection here. The constant k III does not vary appreciably across runs and appears far more constrained than the range reported by Huai et al (2012). The measured and modeled (from the assumed l m and measured u) turbulent stress profiles used in the estimation of D t are also shown in figure B2 for completeness.

Numerical scheme
The numerical runs are conducted in 2D (x-z plane) using the simplified advection-dispersion equation discretized with a central finite difference scheme given as where i = 1, 2 . . . m and j = 1, 2 . . . n represent the node numbers in the x and z directions, respectively. This equation can be rearranged into (m × n) algebraic equation groups for each node given by, where C = {C 1 , C 2 , . . . C mn } represent all the nodes from equation (C6) after rearrangement, and F = {f 1 , f 2 , . . . f mn } represents all the (m × n) equations described by equation (C6). The marching in C can be achieved by using the Newton-Raphson method, where p is the iteration number, C p+1 is the updated solution vector and J = dF/dC is a five-banded Jacobin matrix containing all the coefficients when Figure E4. Variation of removal efficiency e as a function of the dimensionless numbers from the numerical model runs. Figure F5. The coefficients of determination R 2 between c1, c2 and εg . The dashed lines denotes the linear fitting for the 8 FTWs dimensions considered in the model runs.
solving C. For boundary conditions, a unitary concentration is specified at the inlet and a zerohorizontal flux gradient at the outlet so that: The outlet boundary condition is selected as a zero-longitudinal mass flux gradient to ensure that the exit contaminant mass flux from a FTW enters the treated water collection system without modifications. In the vertical direction, diffusive fluxes (the only ones that are solved because mean subsidence is zero) are set to 0 at the bottom of the channel and at the free water surface respectively. A 1-D description of u(z) necessitates a mean advective vertical velocity to be zero assuming water is incompressible and no-slip at the wall-bottom of FTW. These boundary conditions ensure that the FTW outlet is the only exit of contaminants. More complex mechanisms such as contaminant loss to the FTW floor, evaporation or other water losses from the FTW, dissolved gas escape to the atmosphere (due to ebullition, diffusion, or even through the vegetation elements) and so on can be accommodated but are not carried out here to maintain focus on the most generic processes affecting e common to all FTWs.

Appendix D. Longitudinal concentration distribution
To illustrate the determination of e, the longitudinal variation of C H (x) across each vertical section inside the FTW is shown in figure D3. This figure shows the relative vertically integrated concentration (C H (x)/C o ) over the longitudinal (streamwise) direction, which is used in the determination of e. Figure D3 illustrates that the r = 1 representation is more efficient in removing contaminants when compared with the r = 2 (as expected). However, the difference becomes more obvious when small removal kinetic rates (k v < 10 −3 ) are used in the model evaluation.
When the removal kinetic rate becomes high, blockage remains but becomes less significant in terms of influencing e.

Appendix E. Removal efficiency and dimensionless parameters
To formulate a simplified model for removal efficiency, e is examined via four dimensionless parameters: a Da L , a Pe L , ε = H/L and ε g = h/H. The removal efficiency and the dimensionless parameters are plotted in figure E4. Figure E4 shows that Da L has the strongest correlation with e among all the dimensionless parameters and that the r = 1 model is more efficient than the r = 2 in removing contaminants.

Appendix F. Model coefficients: c 1 and c 2
By comparing the FTWs model in equation (10) with the simplified plug flow model, it is reasonable to argue that coefficients c 1 and c 2 should reflect the difference resulting from physical configurations and canopy properties between the two systems. Hence, the correlation of c 1 and c 2 versus ε g = h/H is now examined (since the canopy density a is constant for all the four different ε g cases). The comparison is shown in figure F5. Figure F5 illustrates the correlation between the coefficients c 1 and c 2 and ε g is significant. The coefficient of determinations (R 2 ) are calculated between the dashed lines and optimised values of c 1 and c 2 , which demonstrates that c 1 is significantly correlated (R 2 = 0.99 and 1.00 for r = 1 and r = 2, respectively) with h/H. The c 2 is also closely correlated to h/H (R 2 = 0.64 and 0.73). However, the relation is not as strong as that between c 1 and ε g . Finally, the optimized values for c 1 and c 2 used in the model comparison are summarised in table F2.