1. Introduction
The importance of axisymmetric wake flows lies in their pivotal role in optimising the efficiency and environmental impact of various engineering systems, ranging from wind energy to aerospace design and pollution control. The study of turbulence in axisymmetric wake flows at high Reynolds numbers is a fundamental fluid-mechanics problem treated in many turbulence textbooks (e.g. Pope Reference Pope2000). It is also practically significant because it helps researchers predict and mitigate issues related to drag, flow instability, heat exchange and energy efficiency. For instance, turbulence in wind turbine wakes serves as a double-edged sword. It aids the wake recovery, thus providing more energy for downwind turbines in a wind farm, while simultaneously amplifying the fatigue loads experienced by those very turbines (Stevens & Meneveau Reference Stevens and Meneveau2017).
The need for fast-running models to predict wake flows such as those of wind turbines is essential, as high-fidelity computational fluid dynamics (CFD) modelling or experiments are still too expensive/time consuming and are, therefore, impractical for the optimisation and real-time control of these flows (Meneveau Reference Meneveau2019). Typically, CFD modelling entails choosing a computational discretisation technique, creating often complex numerical grids, selecting a suitable turbulence model and conducting post-processing on the data (Meneveau Reference Meneveau2019). This demands considerable manpower and expertise in fluid mechanics (and turbulence modelling), which typically is not always readily available in industries such as wind energy due to its multidisciplinary nature and diverse workforce. Hence, engineering wake models are often preferred over high-fidelity CFD models in industrial applications. While a substantial amount of research has been dedicated to developing fast-running models to predict the mean velocity distribution in wind turbine wakes (see Porté-Agel, Bastankhah & Shamsoddin (Reference Porté-Agel, Bastankhah and Shamsoddin2020), and references therein), estimation of turbulence in these flows predominantly relies on purely empirical methods (e.g. Crespo & Hernandez Reference Crespo and Hernandez1996; Ishihara & Qian Reference Ishihara and Qian2018) owing to the intricate nature of turbulence modelling. Existing analytical engineering wind-farm models (e.g. Niayifar & Porté-Agel Reference Niayifar and Porté-Agel2016; Zong & Porté-Agel Reference Zong and Porté-Agel2020a; Bastankhah et al. Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021; Lanzilao & Meyers Reference Lanzilao and Meyers2022) typically determine the wake recovery rate based on the incoming turbulence level for each turbine. The turbulent kinetic energy (TKE) generated by upwind turbines, predicted based on empirical models such as Crespo & Hernandez (Reference Crespo and Hernandez1996), then serves as an initial condition for predicting the wake recovery rate of downwind turbines. This emphasises the importance of advancing physics-based methods to predict the TKE distribution. To the best of our knowledge, this work is the first attempt to develop a physics-based engineering TKE model for wake flows. The model derivation is elaborated in § 2, the experimental set-up for two cases used to validate model predictions is described in § 3. Results are then discussed in § 4, before a summary is provided in § 5.
2. Mathematical model development
The axisymmetric wake is described using a cylindrical coordinate system with 
$(x,r,\theta )$, where 
$x$ is the streamwise distance from the object causing the turbulent wake, 
$r$ is the radial distance from the centre of the wake and 
$\theta$ is the azimuthal angle. Mean and fluctuating velocity components in the 
$(x,r,\theta )$ coordinate system are indicated by 
$(U,V,W)$ and 
$(u,v,w)$, respectively. Due to the assumption of axisymmetry, 
$\partial /\partial \theta =0$. In the following, 
$\langle \rangle$ denotes time averaging. The TKE denoted by 
$k$ is defined by 
$k=0.5\langle q^2\rangle$, where 
$q^2=u^2+v^2+w^2$. The steady-state TKE transport equation, neglecting pressure–velocity covariance and viscous diffusion but including swirl, reads as (Shiri Reference Shiri2010)
\begin{align} &\underbrace{U\frac{\partial k}{\partial x}+\left[V\frac{\partial k}{\partial r}\right]}_{\text{advection}}+ \underbrace{\left[\frac{\partial\langle uq\rangle}{2\partial x}\right]+\frac{1}{2r}\frac{\partial \left(r\langle vq \rangle\right)}{\partial r}}_{\text{diffusion}}+ \underbrace{\varepsilon\vphantom{\left(\frac{2{\rm \pi}}{365}\right)}}_{\text{dissipation}}\nonumber\\ &\quad = \underbrace{\left[\langle vw\rangle\frac{W}{r}\right]-\left[\langle w^2\rangle \frac{V}{r}\right]}_{\text{production}}\nonumber\\ &\qquad \underbrace{-\left[\langle v^2\rangle \frac{\partial V}{\partial r}\right]-\left[\langle u^2\rangle \frac{\partial U}{\partial x}\right]-\left[\langle uv\rangle \frac{\partial V}{\partial x}\right]- \left[\langle uw\rangle \frac{\partial W}{\partial x}\right]-\left[\langle vw\rangle \frac{\partial W}{\partial r}\right] - \langle uv\rangle \frac{\partial U}{\partial r}}_{\text{production}}. \end{align}The following simplifications and assumptions are made to be able to mathematically solve (2.1). First, we model the diffusion (i.e. transport) terms based on the gradient-diffusion hypothesis given by (Pope Reference Pope2000)
\begin{equation} \frac{\partial \langle uq\rangle}{2\partial x}+\frac{\partial (r\langle vq\rangle)}{2r\partial r}={-}\frac{\partial}{\partial x}\left(\frac{\nu_T}{\sigma_k}\frac{\partial k}{\partial x}\right)-\frac{1}{r}\frac{\partial}{\partial r}\left( \frac{\nu_T}{\sigma_k}r\frac{\partial k}{\partial r} \right), \end{equation}
where 
$\nu _T$ is the turbulent viscosity and 
$\sigma _k$ is the turbulent Prandtl number. The value of 
$\sigma _k$ generally depends on atmospheric stability (see Li (Reference Li2019) and Basu & Holtslag (Reference Basu and Holtslag2021), among others), with experimental results of 
$\sigma _k \in [0.7,0.92]$ for a neutral atmospheric boundary-layer flow (Businger et al. Reference Businger, Wyngaard, Izumi and Bradley1971; Kays Reference Kays1994). This work assumes 
$\sigma _k=1$ based on the Reynolds analogy, commonly used for turbulent flows (Pope Reference Pope2000). The validity of the gradient-diffusion hypothesis is further examined in § 4. The terms in (2.1) shown in square brackets are normally small compared with other terms for axisymmetric wake flows and can thus be neglected (Uberoi & Freymuth Reference Uberoi and Freymuth1970). Neglecting these terms will be also supported by the budget analysis performed in § 4.
Prior studies (e.g. Wygnanski & Fiedler Reference Wygnanski and Fiedler1970; Hussein, Capp & George Reference Hussein, Capp and George1994) showed that 
$\nu _T$ is fairly uniform at the centre of axisymmetric wakes, but it decays at wake edges. While cross-stream variations of turbulent viscosity can be determined based on velocity profiles (Basset et al. Reference Basset, Viggiano, Barois, Gibert, Mordant, Cal, Volk and Bourgoin2022), for simplicity, the common assumption of 
$\nu _T\approx \nu _T(x)$ is used herein (Pope Reference Pope2000). Moreover, the dominant advection term (i.e. 
$U\partial k/\partial x$) is linearised by replacing 
$U$ with 
$U_0$, where the subscript 
$_0$ denotes the inflow. This approximation improves with distance from the origin, as the velocity deficit 
$\Delta U$ decreases. The TKE dissipation rate is written as
\begin{equation} \varepsilon=C_\varepsilon k^{3/2}/l_m, \end{equation}
where 
$l_m=l_m(x)$ is the mixing length. The normalised energy dissipation rate, 
$C_\varepsilon$, is traditionally assumed to be constant for high Reynolds number flows, but there has been a great deal of evidence in recent years showing that it may not be constant (see the review of Vassilicos (Reference Vassilicos2015), and references therein). Therefore, we assume 
$C_\varepsilon =C_\varepsilon (x)$. Moreover, the turbulent viscosity is commonly modelled by
\begin{equation} \nu_T=ck^{1/2}l_m, \end{equation}
where 
$c$ is a constant (Pope Reference Pope2000). Hence, from (2.3) and (2.4), the TKE dissipation rate 
$\varepsilon$ can be expressed as
\begin{equation} \varepsilon=\frac{\nu_T k}{\varPsi(x)}, \quad \text{where} \ \varPsi(x)=\frac{cl_m^2(x)}{C_\varepsilon(x)}. \end{equation}Finally, the Boussinesq hypothesis is used to model the Reynolds shear stress, where
\begin{equation} \langle uv\rangle ={-}\nu_T \left( \frac{\partial U}{\partial r}+\frac{\partial V}{\partial x}\right). \end{equation}
Note that 
$\partial V/\partial x \ll \partial U/\partial r$, especially in the far-wake region, and thus it can be neglected. Given the hypotheses discussed above, (2.1) can be reduced to
\begin{equation} \frac{U_0}{\nu_t(x)}\frac{\partial k_w(x,r)}{\partial x}-\frac{1}{r}\frac{\partial}{\partial r}\left( r\frac{\partial k_w(x,r)}{\partial r} \right) + \frac{1}{\varPsi(x)}k_w(x,r) \approx \left(\frac{\partial U(x,r)}{\partial r}\right)^2. \end{equation}
Note that, in (2.7), the total TKE, 
$k$, is substituted with the wake-generated TKE, 
$k_w$, where 
$k_w=k-k_0$. This is only possible assuming that the spatial variations of 
$U_0$ and 
$k_0$ are negligible compared with variations caused by the wake. Moreover, the TKE dissipation rate in the background flow is assumed to be considerably smaller than the one in the wake. We therefore neglect the terms including 
$k_0$ on the left-hand side of (2.7), and 
$\partial U/ \partial r$ is only due to the wake-generated shear (i.e. 
$\partial U/\partial r=\partial U_w/\partial r$). Note that the velocity profile 
$U(x,r)$ is assumed to be any smooth function with an arbitrary shape.
The solution of the above equation is sought for the domain of 
$x\geqslant x_0$ and 
$r\geqslant 0$, where 
$x_0$ is the virtual origin. The initial and boundary conditions are defined as 
$k_w(x_0,r)=0$, 
$\partial k_w(x,0)/\partial r=0$ and 
$k_w(x,\infty )\rightarrow 0$ as 
$r\rightarrow \infty$, respectively. Due to the shear in the wake flow, the turbulence level normally starts increasing right from the origin of the wake, so we assume that 
$x_0=0$ is a good approximation. However, 
$x_0$ has been shown to be both positive and negative (Neunaber, Hölling & Obligado Reference Neunaber, Hölling and Obligado2022a), and, if included as a variable parameter, may improve wake model predictions (Neunaber, Hölling & Obligado Reference Neunaber, Hölling and Obligado2024). Thus, 
$x_0$ has an arbitrary value in our model derivation for more flexibility. The solution involves two positive, monotonic functions of 
$x$ and a dummy variable 
$X$ (such that 
$x_0\leqslant X\leqslant x$), namely
\begin{equation} \phi(X,x)=\frac{1}{U_0}\int_{\xi=X}^x\nu_t(\xi) \,\mathrm{d}\xi,\quad \psi(X,x)=\frac{1}{U_0}\int_{\xi=X}^x\frac{\nu_t(\xi)}{\varPsi(\xi)} \,\mathrm{d}\xi, \end{equation}
where 
$\xi$ is a dummy variable. The exact solution of (2.7), achieved using the Green's function method, is
\begin{align} k_w(x,r)&=\int\nolimits_{X=x_0}^x\int\nolimits_{\rho=0}^\infty\frac{\nu_t(X)}{2U_0\phi(X,x)}\nonumber\\ &\quad \times\exp\left\{-\frac{r^2+\rho^2}{4\phi(X,x)}-\psi(X,x)\right\}{\rm I}_0\left(\frac{r\rho}{2\phi(X,x)}\right)\left(\frac{\partial U(X,\rho)}{\partial \rho}\right)^2\rho\,\mathrm{d}\rho\,\mathrm{d}\kern0.7pt X, \end{align}
where 
${\rm I}_0$ is the modified Bessel function of the first kind, and 
$\rho$ is a dummy variable. See Appendix A for more information on the derivation of (2.9). To compute the integrand in (2.9), the wake velocity profile 
$U(x,r)$ needs to be known. This model can be used with any axisymmetric velocity profile, however, since Gaussian-type models are often used to represent the mean flow properties in wakes, these are presented here. A wake with a Gaussian velocity profile is given by (Tennekes & Lumley Reference Tennekes and Lumley1972; Vermeulen Reference Vermeulen1980; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2014)
\begin{equation} U(x,r)=U_0\left[1-C(x)\exp\left(-\frac{r^2}{2\sigma(x)^2}\right)\right], \end{equation}
where 
$\sigma (x)$ is the characteristic wake width at 
$x$, and 
$C(x)$ is the maximum normalised velocity deficit at each 
$x$. A double-Gaussian velocity profile is given by (Schreiber, Balbaa & Bottasso Reference Schreiber, Balbaa and Bottasso2020)
\begin{equation} U(x,r)=U_0 \left[ 1- \frac{1}{2}C(x)\left(\exp\left(-\frac{(r-r_0)^2}{2 \sigma(x)^2}\right)+\textrm{exp}\left(-\frac{(r+r_0)^2}{2\sigma(x)^2}\right) \right) \right], \end{equation}
where 
$r_0$ is the radial position of the double-Gaussian extrema. For a wake with a single-Gaussian velocity-deficit profile as in (2.10), one can simplify (2.9) to
\begin{align} k_w(x,r)&=\int_{X=x_0}^x\frac{U_0\nu_t(X)C^2(x)}{(\sigma^2(X)+4\phi(X,x))^3} \left\{\sigma^2(X)r^2+4\phi(X,x)\left(\sigma^2(X)+4\phi(X,x)\right)\right\}\nonumber\\ &\quad \times\exp\left\{-\frac{r^2}{\sigma^2(X)+4\phi(X,x)}-\psi(X,x)\right\}\mathrm{d}\kern0.7pt X. \end{align}For other wake flow profiles such as double Gaussian, the double integral in (2.9) cannot be reduced to a single integral, so both integrals need to be computed numerically. The reader is referred to Appendix B for a discussion on numerical integration of (2.9) for an arbitrary velocity profile.
3. Experimental set-up
Experiments were conducted in the EnFlo wind tunnel at the University of Surrey with a test section of 
$20\ {\rm m} \times 3.5\ {\rm m}\times 1.5\ {\rm m}$ (
${\rm length} \times {\rm width} \times {\rm depth}$). Tests were run at the free-stream speed 
$1.5\ {\rm m\ s}^{-1}$, as measured by an ultrasonic anemometer mounted at the tunnel inlet (see figure 1). Two different canonical flows are considered in this work: (i) a porous disk in a uniform shear-free flow (i.e. a free-stream (FS) case), and (ii) a turbine model in shear flow over a rough wall (i.e. a boundary-layer (BL) case). These are depicted in figures 1(a) and 1(b), respectively. For the FS case, a porous disk (of diameter 
$D=416$ mm) was manufactured out of a metallic grid arranged in an orthogonal coordinate system, with wires of diameter 1.2 mm and a square grid 
$3.8\times 3.8\ {\rm mm}^2$ in size. The disk was mounted in the centre of the tunnel (both vertically and laterally) via a roof-mounted sting; its estimated thrust coefficient, 
$C_T$, is 0.65. The mean incoming flow for the FS case was verified to be uniform within 
$0.39\,\%$ and characterised by a turbulence intensity of 
$TI=1/U_0\sqrt {1/3(\langle u^2_0\rangle +\langle v^2_0\rangle +\langle w^2_0\rangle })=1.98\,\%$.
Figure 1. Sketch, and overall dimensions, of the wind-tunnel set-up: (a) disk in the free-stream case, (b) wind turbine (WT) in the boundary-layer case and (c) close-up of the three-dimensional laser Doppler anemometry (LDA) set-up (not to scale). Here, 
$\tilde {U}, \tilde {V}, \tilde {W}$ indicate the mean velocities in the Cartesian coordinate system.
For the BL case, a typical offshore BL was achieved by a combination of 13 truncated triangular spires located at the tunnel inlet with roughness elements arranged in a staggered layout covering the entire tunnel floor. The spires (flat plates) have the following dimensions: 60 mm (base), 4 mm (apex), are 600 mm tall and spaced 266 mm in the tunnel's span. The roughness elements were arranged in a staggered layout covering the tunnel floor. The values of 
$(u_*/U_{0})^2$ and roughness length 
$z_0$ are 
$2.2\times 10^{-3}$ and 0.18 mm, respectively, where 
$u_*$ is the friction velocity. The boundary-layer set-up is depicted in figure 1(b), with further information provided in Placidi, Hancock & Hayden (Reference Placidi, Hancock and Hayden2023). The flow mean uniformity across the spanwise direction of the BL at hub height 
$z_h$ is within 
$1.16\,\%$. The wind turbine used in this case is a 
$1:300$ scaled rotating model of a 5 MW offshore wind turbine with a rotor diameter 
$D$ of 
$416$ mm (matching the porous disk case) and a hub height 
$z_h$ of 300 mm. The blades are tapered and twisted flat plates to account for the much lower operating laboratory Reynolds number compared with the full-scale counterpart. See Hancock & Pascheke (Reference Hancock and Pascheke2014) and Placidi et al. (Reference Placidi, Hancock and Hayden2023) for more information on turbine design and characteristics. The model tip-speed ratio 
$\lambda$ was kept at 
$6\pm 1.5\,\%$, which resulted in an estimated 
$C_T$ of 
$0.48$, originally based on wake measurements in uniform flow (Hancock & Pascheke Reference Hancock and Pascheke2014), but later verified with floating-element force balance measurements. The spanwise-averaged hub-height incoming velocity 
$U_0$ is 
$1.42\ {\rm m\ s}^{-1}$, and the incoming turbulence intensity, 
$TI$, is 4.8 %. The turbine was positioned 10 m downstream of the tunnel inlet, as shown in figure 1(b), where the BL had time to fully adjust to the surface conditions and reached a fully developed state, and quantities are slowly varying in 
$x$. A summary of the main experimental parameters for both cases is presented in table 1. Here, the integral time scale is evaluated by integration of the autocorrelation coefficient of the velocity fluctuations until its first zero crossing, with a robust procedure similar to that described in Smith et al. (Reference Smith, Neal, Feero and Richards2018), which involves ensemble averaging the time scale over several independent realisations of the original signal. As in Gambuzza & Ganapathisubramani (Reference Gambuzza and Ganapathisubramani2021), we used signals with a duration of 200 times the initially estimated integral time scale. Then, the non-dimensional spanwise-averaged integral length scale at hub height (
$\varLambda /D$ in table 1) is obtained by applying Taylor's hypothesis of ‘frozen turbulence’.
Table 1. Summary of the experimental conditions in the reference cases (i.e. with no turbine/disk). Boundary-layer characteristics are derived at 
$x/D=2$, although slowly varying in 
$x$. Here, 
$u_*$ is the friction velocity, 
$z_0$ and 
$D$ are the roughness length and the turbine/disk diameter in mm, 
$Re_\delta$ is the BL thickness Reynolds number, 
$\varLambda$ the integral length scale and 
$TI$ is the turbulence intensity.
For both cases, three-component velocity measurements were acquired with a LDA (Dantec Dynamics, Denmark) at a minimum frequency of 200 Hz for 120 s for each measurement point, which balanced competing requirements between data statistical convergence, reasonable running time, seeding density and temporal resolution required to resolve the turbulence scales of interest, as further discussed in Placidi et al. (Reference Placidi, Hancock and Hayden2023). The three components are acquired independently, but can be synchronised by interpolation, when required by the analysis. Standard errors are within 
$\pm 0.5\,\%$ and 
$\pm 5\,\%$ for the mean and second-order quantities (95 % confidence level). Three laser beams emanating from two laser probes (of a focal length of 
$300$ mm) were used in conjunction to measure the velocity components as shown in figure 1(b). One probe measured streamwise and lateral components, while the other measured the vertical component independently. A 
$45^{\circ }$ mirror is used to focus the beam measuring the vertical component onto the same measurement volume of the other two beams, allowing for simultaneous measurements of all three velocity components. Both probes were mounted vertically in the tunnel (hence the need for the mirror) and embedded into an aerodynamic shroud to minimise flow interference. This set-up helps minimise the intrusiveness of the measurement system while circumventing the error propagation originated by the 3-D transformation matrix and accurate determination of the position/separation between the beams. Measurements, for all cases, were collected at different streamwise locations (
$2\leqslant x/D\le 15$) both with and without the turbine/disk model to isolate the wake-added quantities from their counterpart in the background flow.
Before any results are presented, we discuss the relevance of the BL case to the model developed for axisymmetric wake flows in § 2. The TKE transport equation (in § 2) is simplified by assuming wake flow axisymmetry and inflow homogeneity. While these assumptions are valid for the porous disk in a uniform flow, they do not hold in the BL case. Figure 2(a) shows vertical profiles of normalised streamwise velocity (
$U/U_0$) and TKE (
$k/U_0^2$) at two downwind locations for the BL case, where 
$z$ is the height from the ground in the Cartesian coordinate system. Inflow conditions are also reported for comparison (dashed lines). Figure 2(a) shows that, due to the mean shear in the incoming BL, neither the inflow velocity nor the inflow TKE is uniform in the vertical direction. Looking at 
$x=5D$, it is evident that the wake increases the flow shear in the upper half of the wake while decreasing it in the lower half. This generates/suppresses turbulence in above/below hub height, as shown in figure 2(a). The assumptions made in § 2 are therefore clearly violated in the vertical direction, and the developed model is not expected to provide satisfactory predictions. In the lateral direction, 
$y$, however, the incoming flow is approximately uniform, and the TKE lateral profiles appear to be more symmetrical, as seen in figure 2(b). Therefore, despite the fact that the wake is not axisymmetric in this case, it is still of interest to examine whether the developed model can be employed to estimate lateral TKE profiles at the turbine's hub height. It is also worth noting that the slight lateral wake deflection towards the negative y direction, as seen in figure 2(b), is due to the interaction of the wake swirling in the anticlockwise direction (seen from upstream) with the incoming flow shear (Fleming et al. Reference Fleming, Gebraad, Lee, van Wingerden, Johnson, Churchfield, Michalakes, Spalart and Moriarty2014). Since equations in § 2 are written in cylindrical coordinates, hereafter, quantities for the BL case are averaged on both sides of the wake for presentation purposes (i.e. 
$f(r)=0.5(f(y)+f(-y))$ for 
$r=y$).
Figure 2. (a) Vertical and (b) lateral profiles of normalised streamwise velocity (
$U/U_0$) and normalised TKE (
$k/U_0^2$) for the BL case. Horizontal lines represent the hub and the tip positions.
4. Results and discussions
To predict TKE profiles based on (2.9), the wake velocity profile 
$U(x,r)$ needs to be first known either from experiments or engineering wake models. Figure 3 shows radial profiles of the normalised streamwise velocity deficit, 
$\Delta U/U_0=(U-U_0)/U_0$, at different streamwise locations for both cases and compares them with two customary profiles: Gaussian (2.10) and double Gaussian (2.11). The figure shows that the velocity deficit in the BL case is consistently smaller than the other case. The observed difference is expected to be caused by the different inflow conditions and thrust coefficients in the two cases. A higher level of inflow turbulence and a lower value of thrust coefficient in the BL case leads to a less pronounced wake. To predict wake velocity profiles, the double-Gaussian profile (2.11) is used hereafter as it generally better represents our experimental data for 
$x/D \leqslant 6$ where the wake profiles present a double peak, although some discrepancies are noticeable, especially for the FS case. Further work can be done to see if the double-Gaussian model parameters can be better optimised to improve its predictions. It is important to note that inaccuracies in velocity prediction in the near wake manifest as an error in TKE predictions for the far wake. This highlights the importance of accurate velocity predictions in the near wake. Depending on the actual shape of the near-wake profiles, one may also apply other wake models such as super-Gaussian (Shapiro et al. Reference Shapiro, Starke, Meneveau and Gayme2019; Blondel & Cathelain Reference Blondel and Cathelain2020). This is possible due to the versatility of the new model developed for a generic profile of 
$U(x,r)$.
Figure 3. Radial profiles of normalised velocity deficit (
$\Delta U/U_0$) at different 
$x/D$ based on experiments, Gaussian and double Gaussian (referred to as D-Gaussian) wake models. Solid lines show the BL case and dashed lines show the FS case.
Next, the radial profiles of the normalised azimuthal velocity 
$W/U_0$ at different streamwise positions are shown in figure 4(a) for both cases. The azimuthal velocity in the FS case is negligible as the porous disk does not generate any swirl motion. In the BL case, while rotating blades of the turbine induce swirl motion in the wake, the value of wake swirl decays rapidly. Figure 4(b) shows that the maximum azimuthal velocity in the wake is considerably smaller than the maximum velocity deficit, especially in the far-wake region. This is consistent with the TKE budget analysis, discussed later, which will show that all terms in the TKE budget (2.1) that include azimuthal velocity have negligible values at 
$x=5D$. Nonetheless, it is still important to bear in mind that, despite the rapid decay of swirl and seemingly its small contribution in the far wake, the wake swirl in the near wake may still have a non-negligible impact on the TKE in the far wake.
Figure 4. (a) Lateral profiles of normalised azimuthal velocity (
$W/U_0$) at hub height. (b) Ratio of the maximum 
$W$ to the maximum 
$\Delta U$ at different 
$x/D$.
Figure 5 shows radial profiles of the wake-added turbulent velocity fluctuations (
$\langle u^2 \rangle _w$, 
$\langle v^2 \rangle _w$, and 
$\langle w^2 \rangle _w$) and the wake-added TKE (
$k_w$) at different streamwise locations for both cases. It is clear how the wake edge (i.e. 
$r\approx 0.5 D$) where the shear production is maximum corresponds to the maximum level of turbulence. At the wake edge, 
$\langle u^2 \rangle _w$ is dominant and significantly larger than 
$\langle v^2 \rangle _w$ and 
$\langle w^2 \rangle _w$ in both cases. This is especially the case in the near wake of the BL case. However, the other two components (
$\langle v^2 \rangle _w$, and 
$\langle w^2 \rangle _w$) gradually become larger and more comparable to 
$\langle u^2 \rangle _w$ in the far wake. Spalart (Reference Spalart1988) showed that, for BL flows, the TKE diffusion due to fluctuating pressure is mainly small, however, pressure has a significant role in redistributing the energy by extracting it from the streamwise component and transferring it to the other two components (Pope Reference Pope2000). This seems to be the case here too. It is also interesting to note that, at the wake centre (
$r/D<0.25$) of the BL case, 
$\langle v^2 \rangle _w$ and 
$\langle w^2 \rangle _w$ are bigger than 
$\langle u^2 \rangle _w$ for 
$x>3D$.
Figure 5. Radial profiles of normalised wake-added turbulent quantities at different 
$x/D$. Solid lines show the BL case and dashed lines show the FS case.
Cross-stream turbulent fluctuations are especially critical in explaining unsteady oscillations of wakes termed as wake meandering (Larsen et al. Reference Larsen, Madsen, Thomsen and Larsen2008), which has major impacts on flow mixing and wake expansion. However, prior experiments have often quantified the TKE distribution only based on streamwise fluctuations. Moreover, steady-state engineering wake models have often used the turbulence intensity defined only based on streamwise fluctuations to estimate the wake expansion. This highlights the importance of the current experimental work in capturing the total TKE.
Next, we determine three important turbulent quantities using the experimental data: (i) turbulent viscosity (
$\nu _t$), (ii) mixing length (
$l_m$) and (iii) the normalised TKE dissipation rate (
$C_\varepsilon$). The last two are needed to compute the parameter 
$\varPsi =cl_m^2/C_\varepsilon$ in the simplified TKE transport equation (2.7). Both turbulent viscosity and mixing length are estimated based on the method described in Bai, Meneveau & Katz (Reference Bai, Meneveau and Katz2012) and later implemented in Rockel et al. (Reference Rockel, Peinke, Hölling and Cal2016) and Scott et al. (Reference Scott, Martínez-Tossas, Bossuyt, Hamilton and Cal2023). At each streamwise location, 
$\nu _t$ is the slope of the linear curve fitted to the variation of 
$-\langle uv\rangle$ with respect to 
$\partial U/\partial r$, and 
$l_m^2$ is the slope of the linear curve fitted to the variation of 
$-\langle uv\rangle$ with respect to 
$|\partial U/\partial r|\partial U/\partial r$.
Figure 6(a) shows that, in both cases as expected, the normalised mixing length increases with 
$x/D$, which indicates an increased characteristic length scale as the wake flow evolves (Iungo et al. Reference Iungo, Viola, Ciri, Rotea and Leonardi2015). The trend is almost identical for the FS case, but with a lower value at each 
$x/D$, which would be expected because of the lower value of inflow integral length scale in this case (see table 1). The normalised turbulent viscosity, shown in figure 6(b), also linearly increases with 
$x$ and seems to approach an almost constant value in the far wake. The same trend is shown in experiments by Zong & Porté-Agel (Reference Zong and Porté-Agel2020b), displayed in green. In all three cases shown in figure 6(b), the turbulent viscosity reaches a plateau at approximately 
$6D-8D$ downstream. Experiments and numerical simulations by Scott et al. (Reference Scott, Martínez-Tossas, Bossuyt, Hamilton and Cal2023) (not shown here) also reported a fairly similar behaviour, and found that, further downstream (e.g. 
$x/D\gg 15$), the turbulent viscosity should decrease as the wake recovers. This behaviour of 
$\nu _t$ can be explained with the hypothesis that 
$\nu _t \approx l_m^2|\partial U/\partial r|$ (Bai et al. Reference Bai, Meneveau and Katz2012). In the near wake, the increase in the turbulent viscosity is mainly due to the increase in the mixing length. However, in the far wake, the growth in the mixing length is balanced by the decreasing velocity gradients across the wake. It is also noteworthy that, despite having higher 
$C_T$, the turbulent viscosity is smaller in the FS case compared with the BL case, which can be attributed to a lower level of inflow turbulence in the former case. Furthermore, the turbulent viscosity in Zong & Porté-Agel (Reference Zong and Porté-Agel2020b) is higher than the BL case although both have a fairly similar inflow turbulence level. The discrepancy between these two cases can be explained by the different thrust coefficients of the turbine models. In conclusion, figure 6(b) shows that, while the variation of 
$\nu _t$ with 
$x$ follows a fairly similar pattern in all cases, its value is sensitive to both inflow turbulence and thrust force.
Figure 6. Evolution of normalised (a) mixing length (
$l_m/D$), (b) turbulent viscosity (
$\nu _t/U_0D$) including data reported by Zong & Porté-Agel (Reference Zong and Porté-Agel2020b), (c) energy dissipation rate (
$C_\varepsilon$) and (d) 
$l_m^2/C_\varepsilon D^2$, with 
$x/D$. The black circles are from the BL case and the blue circles are from the FS case.
Next, we compute the TKE dissipation rate (
$\varepsilon$, hereafter dissipation in short). Estimating the dissipation from experimental data is a troublesome task (Wang et al. Reference Wang, Yang, Wu, Ma, Peng, Liu and Wang2021). Two common approaches are used here. In the first method, thanks to our comprehensive experimental data set, we estimate the dissipation indirectly by performing a TKE budget analysis, i.e. evaluating all the terms in (2.1), so that by definition its residual is the dissipation (Hearst & Lavoie Reference Hearst and Lavoie2014). As an example, the TKE budget is provided in figure 7 for 
$x=5D$. Figure 7 demonstrates that all the terms in square brackets in (2.1) can, indeed, be considered negligible for both cases. The diffusion terms are shown in red, where the real terms are represented as circles and the modelled terms based on the gradient-diffusion hypothesis in (2.2) as lines. It can be seen that, in both cases, the radial diffusion term, which is the dominant diffusion mechanism in the wake, and its modelled counterpart behave similarly. The figure also suggests that the streamwise convection term is in balance with the radial diffusion term at this streamwise location, which is particularly evident in the BL case. The radial diffusion term is positive at the wake edge, acting as a sink of energy transporting energy towards the wake's centre and the outer region where the radial diffusion is negative (i.e. acting as an energy source). This process flattens and widens the TKE profiles as the wake moves downstream, as already seen in figure 5. The dominant production term (dot-dashed black line) is mainly in balance with the dissipation (green dotted line) as the main two dominant terms and, as expected, 
$-\varepsilon$ is negative so that it acts as an energy sink. Both have maximum absolute values in correspondence of the blade tip at the wake edge (
$r/D\approx 0.5$).
Figure 7. The TKE budget for both cases at 
$x/D=5$ (terms normalised by 
$U_0^3/D$). Advection (blue), diffusion (red), production (black), dissipation (green).
Given the time-resolved nature of the experimental data, we can also directly quantify the dissipation based on another common approach that assumes the turbulence at small scales to be homogeneous and isotropic (subscript 
$_{HIT}$). This leads to the estimation of the dissipation based on Taylor's hypothesis of frozen turbulence, as 
$\varepsilon _{HIT}=15\nu \langle (\partial u/\partial t)^2\rangle /U^2$, where 
$\nu$ is the kinematic viscosity and 
$t$ is time (Dairay, Obligado & Vassilicos Reference Dairay, Obligado and Vassilicos2015; Neunaber, Peinke & Obligado Reference Neunaber, Peinke and Obligado2022b). Figure 7 shows there is a satisfactory agreement between the two methods used herein to estimate the dissipation, particularly in the BL case. Although there is a difference in the magnitude between the two dissipation terms at the wake edge (
$r/D\approx 0.5$) for the FS case, overall there is a satisfactory agreement given a high level of uncertainty in estimating the dissipation.
Once the dissipation is estimated, the normalised TKE dissipation rate (
$C_\epsilon$) can be computed based on 
$\epsilon l_m/{K}_w^{3/2}$, where 
$\epsilon$ and 
${K}_w$ are the maximum absolute values of 
$\varepsilon$ and 
$k_w$ at each streamwise location, respectively. Figure 6(c) confirms that 
$C_\varepsilon$ is, indeed, not constant and instead increases with 
$x$. This is in line with previous works that highlighted how 
$C_\varepsilon$ may vary in non-equilibrium flows (Dairay et al. Reference Dairay, Obligado and Vassilicos2015; Vassilicos Reference Vassilicos2015; Obligado, Dairay & Vassilicos Reference Obligado, Dairay and Vassilicos2016). Moreover, in line with other turbulent quantities, the normalised TKE dissipation rate is higher for the BL case. Next, figure 6(d) reports the variation of 
$l_m^2/C_\varepsilon D^2$ with 
$x/D$. Based on the fitted linear curve shown by the black dashed line in the figure, 
$l_m^2/C_\varepsilon D^2\approx 0.0072+0.0032x/D$ for the BL case. As discussed in § 2, 
$\varPsi =cl_m^2/C_\epsilon$, so 
$\varPsi _{BL}/D^2\approx c(0.0072+0.0032x/D)$. Similarly, for the FS case (blue dashed line), 
$\varPsi _{FS}/D^2\approx c(0.0076+0.0039x/D)$, where 
$c=0.46$ seems to provide satisfactory predictions for both. It is interesting to note that this value is comparable to 
$c=0.55$ used in the log layer of boundary-layer flows (Pope Reference Pope2000). Furthermore, despite significant disparities in other turbulent parameters such as 
$l_m$ and 
$C_{\varepsilon }$, the value of 
$\varPsi$ appears fairly similar in both cases. This similarity tempts us to speculate about the existence of a universal relationship for 
$\varPsi$. However, further investigation is required to substantiate this hypothesis.
After assessing both 
$\nu _t(x)$ and 
$\varPsi (x)$, we proceed to compare the predictions of the model outlined in § 2 with the experimental data, depicted in figure 8. Some deviations are apparent in predicting the location of TKE maxima in the near wake for the FS case and also the TKE level at the wake centre in the far wake for both cases. These discrepancies are believed to primarily stem from inaccuracies in the velocity-deficit predictions illustrated in figure 3. Nevertheless, the overall trend reveals that the model reasonably predicts both the magnitude of wake-added TKE (
$k_w$) and its radial distribution for 
$3\leqslant x/D\leqslant 15$ in both cases. Moreover, the model successfully captures the difference in wake-added TKE between the two cases. Figure 8 shows that the wake-added TKE is considerably higher in the FS case compared with the BL case. As previously discussed, our results show that the inflow turbulence and the amount of the thrust force significantly impact all pertinent turbulent quantities. Thus, predicting their overall effects on wake-added TKE is not straightforward, given their counteracting influences on wake behaviour. This emphasises the necessity of employing a physics-based TKE model capable of realistically predicting wake behaviour under varied operating conditions. It is worth highlighting that, while the FS case shows approximately a 
$3\,\%$ lower turbulence intensity compared with the BL case (as indicated in table 1), the difference becomes significant when considering the integral length scale, where there is more than a one order of magnitude difference. This emphasises that relying solely on turbulence intensity may not fully capture the impact of the incoming turbulent flow on the wake evolution (Gambuzza & Ganapathisubramani Reference Gambuzza and Ganapathisubramani2021; Hodgson, Madsen & Andersen Reference Hodgson, Madsen and Andersen2023).
Figure 8. Model predictions of radial profiles of normalised wake-added TKE (
$k_w/U_0^2$) in comparison with the experimental data and the empirical models of Crespo & Hernandez (Reference Crespo and Hernandez1996) and Ishihara & Qian (Reference Ishihara and Qian2018) at different streamwise positions. The BL case is shown in black and the FS case in blue.
Predictions based on the empirical models of Crespo & Hernandez (Reference Crespo and Hernandez1996) and Ishihara & Qian (Reference Ishihara and Qian2018) are also shown in figure 8. Since these two models predict the wake-added turbulence only based on streamwise velocity fluctuations, for a fair comparison, we use the relationship 
$\sqrt {\langle u_w^2\rangle }=\alpha \sqrt {k_w}$, where 
$\alpha$ is a constant (Larsén Reference Larsén2022). Various values for 
$\alpha$ have been suggested in the range 
$\alpha \in [0.82 \ 1.03]$ (e.g. Cleijne Reference Cleijne1992; Crespo & Hernandez Reference Crespo and Hernandez1996; Malki et al. Reference Malki, Masters, Williams and Croft2014). Our experimental data suggest that 
$\alpha \approx 0.93$ and this is used to plot the total wake-added TKE (
$k_w$) predictions for these two empirical models in figure 8. The predictions based on the model of Crespo & Hernandez (Reference Crespo and Hernandez1996) only depends on the streamwise distance from the disk/turbine, while the one proposed by Ishihara & Qian (Reference Ishihara and Qian2018) also predicts the radial distribution of wake-added turbulence. Both empirical models cannot fully capture the impact of changing operating conditions on the wake-added TKE distribution.
5. Summary
A new fast-running model to predict the 3-D TKE distribution in axisymmetric wake flows is presented. Detailed 3-D LDA measurements were conducted for two canonical cases: (i) a porous disk exposed to a uniform FS flow, and (ii) a turbine model under a turbulent BL. While the former configuration generates an axisymmetric wake, the wake flow in the latter case is non-axisymmetric due to inflow shear. In the latter case, our analysis primarily focused on the flow distribution within a horizontal plane at hub-height level. A budget analysis is first performed to identify dominant terms in the TKE transport equation written in a cylindrical coordinate system. The Boussinesq turbulent-viscosity and gradient-diffusion hypotheses were used to simplify production and diffusion terms, respectively. The simplified partial differential equation was then solved using the Green function's method, which led to a solution written in the form of a double integral. Further simplifications were applied to the exact mathematical solution to facilitate the numerical integration. The new model requires numerical integration based on simple methods such as the trapezoid rule, which can be performed with very basic mathematical knowledge. Therefore, in addition to addressing computational costs, the ease of use (in comparison with CFD models) is the main driving factor behind the development of such a model. The developed solution predicts second-order flow statistics (i.e. TKE distribution) from the knowledge of first-order flow statistics (i.e. time-averaged streamwise velocity distribution, 
$u(x,r)$). While the solution was derived for an arbitrary distribution of 
$u(x,r)$, a double-Gaussian profile was assumed herein due to its resemblance to the experimental data, especially in the near wake. To predict the TKE, the model also necessitates the turbulent viscosity 
$\nu _t(x)$, and the parameter 
$\varPsi (x)=cl_m^2/C_\varepsilon$.
The experimental data showed that in both cases the mixing length in the wake flow grows with the streamwise distance from the disk/turbine. The increase in the mixing length initially leads to an increase in turbulent viscosity, but the turbulent viscosity approaches a constant value in the far wake as wake velocity gradients diminish with the wake recovery. Operating conditions (i.e. 
$C_T$ and the inflow turbulence) are found to have major impacts on turbulent quantities such as the mixing length and turbulence viscosity. Moreover, in agreement with non-equilibrium similarity theory (Vassilicos Reference Vassilicos2015), we observed that the normalised TKE dissipation rate (
$C_\varepsilon$) is indeed not constant, but it increases with streamwise distance from the turbine; this results in a linear relationship for 
$\varPsi (x)$. Finally, the new TKE model predictions are compared with the experimental data, demonstrating a satisfactory level of agreement both in the magnitude and the radial shape of the TKE profiles across a wake flow extent of interest to wind-farm developers and operators (
$3< x/D<15$).
This section is concluded by a discussion on model limitations and future research directions. In this work, we relied on the same experimental data to determine the variations of 
$\nu _t(x)$ and 
$\varPsi (x)$ with 
$x$, which were then used to validate model predictions. This validation does not inherently establish the universality of the relationships governing these input parameters. Therefore, we opt not to claim that the same parameter settings are universally applicable across all different scenarios. Our main goal in this study was to demonstrate the robustness of the developed model. Specifically, we aimed to show that when accurately estimating these parameters, our model reliably predicts the distribution of TKE in axisymmetric wake flows (and approximates non-axisymmetric turbine wake flows at the hub-height level). However, we acknowledge the need for further research to establish universal relationships for 
$\nu _t$ and 
$\varPsi$ across a range of relevant parameters, thereby creating a comprehensive framework for TKE engineering modelling.
The developed model uses information on the mean streamwise velocity distribution as an input to estimate TKE generation resulting from flow shear. Integrating this model into existing engineering velocity-deficit models that determine the wake recovery rate based on the incoming turbulence is straightforward. However, this approach does not take into account the two-way coupling effect, where velocity gradients in the wake generate higher turbulence levels, which subsequently influence mean flow distribution through enhanced flow mixing. Investigating how the wake-added TKE may impact flow entrainment and wake recovery presents an interesting area of research (Nygaard et al. Reference Nygaard, Steen, Poulsen and Pedersen2020). Future studies could potentially explore this by using the developed model in an iterative approach that considers the interplay between first-order and second-order statistics. Alternatively, addressing this coupling in the developed TKE model may involve further simplification and decoupling of equations, leveraging assumptions such as the self-similarity of TKE profiles. In addition, more refined models are needed to predict vertical TKE profiles in non-axisymmetric wakes of turbines immersed in BL flows.
Acknowledgements
We would like to thank Dr P. Hayden, who aided the experimental data collection.
Funding
The experimental work was supported via Flex Funding agreement project VENTI within EPSRC-Supergen ORE (EP/S000747/1) coordinated by the University of Plymouth.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The wind-tunnel data utilised in this work are available at the following link: https://doi.org/10.15126/surreydata.901120.
Appendix A. Green's function solution of the TKE transport equation
The derivation of (2.9) is as follows. First, we use an integrating factor to cancel the 
$k_w(x,r)$ term (i.e. the third term on the left-hand side) in (2.7). Set 
$k_w(x,r)=\exp \{-\psi (x_0,x)\}Q(x,r)$, which reduces (2.7) to
\begin{equation} \frac{U_0}{\nu_t(x)}\frac{\partial Q(x,r)}{\partial x}-\frac{1}{r}\frac{\partial}{\partial r}\left( r\frac{\partial Q(x,r)}{\partial r} \right) = \left(\frac{\partial U({x},r)}{\partial r}\right)^2\exp\{\psi(x_0,x)\}. \end{equation}
Now introduce a change of variable 
$T(x)=\phi (x_0,x)$ and a corresponding dummy variable 
$t(X)=\phi (x_0,X)$ as a coordinate on the interval 
$0\leqslant t\leqslant T$. As 
$\nu _t>0$ throughout the region of interest, 
$t$ is a monotonically increasing function of 
$X$; consequently, 
$X$ can be used as a coordinate in the range 
$x_0\leqslant X\leqslant x$.
Let 
$Q(x,r)=q(T(x),r)$, which reduces (A1) to an inhomogeneous axisymmetric heat equation
\begin{equation} \frac{\partial q(T,r)}{\partial T}-\frac{1}{r}\frac{\partial}{\partial r}\left( r\frac{\partial q(T,r)}{\partial r} \right) = \left(\frac{\partial U({x}(T),r)}{\partial r}\right)^2\exp\{\psi(x_0,x(T))\}. \end{equation}
The Green's function for such equations, subject to 
$q(0,r)=0$, 
$\partial q/\partial r(T,0)=0$ and 
$q(T,r)\rightarrow 0$ as 
$r\rightarrow \infty$, is (see Cole et al. Reference Cole, Beck, Haji-Sheikh and Litkouhi2011)
\begin{equation} G(T,r;t,\rho)=\frac{H(T-t)}{4{\rm \pi}(T-t)}\exp\left\{-\frac{r^2+\rho^2}{4(T-t)}\right\}{\rm I}_0\left(\frac{r\rho}{2(T-t)}\right), \end{equation}
where 
$H$ is the Heaviside step function and 
${\rm I}_0$ is the modified Bessel function of the first kind. Integrating over the cylindrical domain yields the solution of (A2)
\begin{align} q(T,r)&=\int_{t=0}^T\int_{\rho=0}^\infty\frac{1}{2(T-t)}\nonumber\\ &\quad\times\exp\left\{-\frac{r^2+\rho^2}{4(T-t)}+\psi(x_0,X(t))\right\}{\rm I}_0 \left(\frac{r\rho}{2(T-t)}\right)\left(\frac{\partial U(X({t}),\rho)}{\partial \rho}\right)^2\rho\,\mathrm{d}\rho\,\mathrm{d}t{{.}} \end{align}
Changing the dummy variable from 
$t$ to 
$X$ yields the solution (2.9) for 
$k_w(x,r)$.
Appendix B. Numerical integration of (2.9)
To ease the numerical integration of (2.9) for an arbitrary wake velocity profile, we apply the changes below to the exact solution:
(i) For large values of
$\rho$, the Bessel function in (2.9) goes to infinity, while the exponential term goes to zero, which may introduce errors in the numerical integration. For $r>0$, it is useful to write the solution (2.9) in the following form:
(B1)The expression in square brackets is a well-behaved function,
\begin{align} k_w(x,r)&=\frac{1}{U_0}\int_{X=x_0}^x\int_{\rho=0}^\infty\nu_t(X)\exp({-\psi(X,x)}) \left\{\frac{\exp\left({-\dfrac{(r-\rho)^2}{4\phi(X,x)}}\right)}{\sqrt{4{\rm \pi}\phi(X,x)}}\right\}\nonumber\\ &\quad\times\left[\sqrt{\frac{{\rm \pi} r\rho}{\phi(X,x)}}\exp\left({-\frac{r\rho}{2\phi(X,x)}}\right){\rm I}_0\left(\frac{r\rho}{2\phi(X,x)}\right)\right] \left(\frac{\partial U(X,\rho)}{\partial \rho}\right)^2\nonumber\\
&\quad\times\sqrt{\frac{\rho}{r}}\,\mathrm{d}\rho\,\mathrm{d}\kern0.7pt X. \end{align}
$f(z)$, where $z=r\rho /\phi (X,x)$, whose graph rises rapidly from $0$ up to around $1.2$ at $z\approx 1.5$, then decreases towards its limiting value of $1$ as $z$ increases further. By combining the power series and asymptotic expansions for the modified Bessel function, one obtains the following good approximation:
(B2)which gives a maximum error of less than
\begin{equation} f(z)\approx\begin{cases} \sqrt{{\rm \pi} z}\exp({-}z/2)\left(1+\dfrac{z^2}{16}+\dfrac{z^4}{1024}\right), & 0\leqslant z \leqslant 4;\\ 1+\dfrac{1}{4z}+\dfrac{9}{32z^2}, & z>4, \end{cases} \end{equation}
$1.3\,\%$.(ii) The integrand in the exact solution (2.9) has a singularity at
$X=x$, because $\phi (x,x)=0$. To avoid this singularity, we compute the integral by restricting $X$ to the interval $[x_0,x-\delta ]$, where $\delta$ is the size of the grid used for numerical integration. This approximation provides satisfactory results with negligible error for small values of $\delta$ (e.g. $\delta \leqslant 0.1D$, where $D$ is the diameter of the object causing the turbulent wake).(iii) In the exact solution (2.9), the upper bound of integration with respect to
$\rho$ is infinity. For the numerical integration, we replace this with a large finite value, namely $3D$. The velocity gradient $\partial U(X,\rho )/\partial \rho$ quickly goes to zero for large values of $\rho$, so this has a negligible effect on final results.
In summary, instead of the exact solution (2.9), one can numerically compute
\begin{align} k_w(x,r)&=\frac{1}{U_0}\int_{X=x_0}^{x-\delta}\int_{\rho=0}^{3D}\nu_t(X)\exp({-\psi(X,x)})\nonumber\\ &\quad \times\left\{\frac{\exp\left({-\dfrac{(r-\rho)^2}{4\phi(X,x)}}\right)}{\sqrt{4{\rm \pi}\phi(X,x)}}\right\} f\left(\frac{r\rho}{\phi(X,x)}\right) \left(\frac{\partial U(X,\rho)}{\partial \rho}\right)^2\, \sqrt{\frac{\rho}{r}}\,\mathrm{d}\rho\,\mathrm{d}\kern0.7pt X, \end{align}
where the function 
$f(z)$ is approximated by (B2). It is worth noting that (B1), and its approximated form (B3), are valid for 
$r>0$. At 
$r=0$, the exact solution (2.9) is simplified to
\begin{align} k_w(x,0)=\int_{X=x_0}^{x-\delta}\int_{\rho=0}^{3D}\frac{\nu_t(X)}{2U_0\phi(X,x)}\exp\left\{- \frac{\rho^2}{4\phi(X,x)}-\psi(X,x)\right\}\left(\frac{\partial U(X,\rho)}{\partial \rho}\right)^2\,\rho\,\mathrm{d}\rho\,\mathrm{d}\kern0.7pt X. \end{align}
Open access
