**2. Numerical Method**

The Virtual Flow Simulator (VFS-Wind) code [23–25] is employed in this work to simulate wakes from the two turbines, which solves the incompressible Navier–Stokes equations in curvilinear coordinates shown as follows:

$$J\frac{\partial \mathcal{U}^j}{\partial \mathcal{J}^j} = 0,\tag{1}$$

$$\frac{1}{f}\frac{\partial \mathcal{U}^i}{\partial t} = \frac{\mathbb{S}^i\_l}{f}\left(-\frac{\partial (\mathcal{U}^j u\_l)}{\partial \mathcal{S}^j} + \frac{\mu}{\rho}\frac{\partial}{\partial \mathcal{S}^l}\left(\frac{\mathcal{S}^{jk}}{f}\frac{\partial u\_l}{\partial \mathcal{S}^k}\right) - \frac{1}{\rho}\frac{\partial}{\partial \mathcal{S}^l}\left(\frac{\mathcal{Z}^j\_l p}{f}\right) - \frac{1}{\rho}\frac{\partial \tau\_{lj}}{\partial \mathcal{S}^l} + \frac{1}{\rho}f\_l\right),\tag{2}$$

where *xi* (*i* = 1, 2, 3) are the Cartesian coordinates, *ξ<sup>j</sup>* (*j* = 1, 2, 3) are the curvilinear coordinates, *U<sup>i</sup>* denotes the contravariant volume flux, *ui* are the velocity vector in Cartesian coordinates, *J* is the Jacobian of the geometric transformation, *ξ<sup>i</sup> <sup>l</sup>* represent the transformation metrics, *<sup>g</sup>jk* represents the contravariant metric tensor, *ρ* is the density, *μ* is the dynamic viscosity, *p* is the pressure, *fl*(*l* = 1, 2, 3) are the body forces from actuator type turbine models, and *τij* is the subgrid scale stress modeled using the dynamic subgrid scale model [26]. The governing equations are discretized in space using a second-order central differencing scheme and in time using a second-order accurate fractional step method. The pressure Poisson equation is solved using an algebraic multigrid acceleration along with GMRES solver. The momentum equation is solved using the matrix-free Newton–Krylov method. The VFS-Wind code has been validated extensively using laboratory and field measurements [21,24] and applied to utility-scale wind turbines [22,27,28].

The actuator surface models for turbine blades and nacelle developed in [25] are employed in this work. In the actuator surface model for turbine blades, the blade is represented with the actuator surface defined by chords at different radial locations. Forces computed using the blade element method are distributed on the actuator surface to represent the effect of blades on the incoming flow. Compared with the actuator line parameterization, the actuator surface can represent better the geometrical effect of the blade in the chordwise direction. A model for the nacelle is also necessary for accurately predicting turbine wakes, which affects the hub vortex and meandering in the far wake [29]. In this work the nacelle is represented using an actuator surface formed by the actual surface of the nacelle. The effects of nacelle on the incoming flow are modeled using distributed forces with the wall-normal component calculated by satisfying the no-flux boundary condition and the wall-tangent component calculated by specifying a friction coefficient.
