## 2 Problem Description  
Figure 1 sketches the computational fluid domain (\(\Omega\)) used in Lagrangian vortex method simulations, which contains a pair of vortical structures shedding from wingtips and rotating in opposite directions. The fluid domain is delimited by the ground plane surface, \(\mathrm{S_{b}}\), and the surface \(\mathrm{S_{so}}\) far from the wall plane, such that \(\mathrm{S=S_{b}\cup S_{so}}\). The boundary conditions for the mathematical formulation must be imposed on \(\mathrm{S}\). The circulation strength is defined by \(\Gamma^{*}=\pm\mathrm{W^{*}}/(\mathrm{pb^{*}U_{a}}^{*})\), where \(\mathrm{W^{*}}\) is the aircraft weight, \(\mathrm{\rho}\) is the fluid density, \(\mathrm{b^{*}}\) is the wingspan, and \(\mathrm{U_{a}}^{*}\) is the aircraft approaching velocity. Here, all dimensional quantities are identified through the symbol \({}^{*}\).  
In Figure 1, the inertial frame origin is located at \(\mathrm{x^{*}}=0.0\) and \(\mathrm{y^{*}}=0.0\). The ground plane length is defined by \(\mathrm{B^{*}}\), and the starting positions of the center of the counter-rotating vortical structures are given by the coordinates \(\mathrm{x^{*}}=\mathrm{B^{*}}/2\)\(\pm\)\(\mathrm{b^{*}}/2\) and \(\mathrm{y^{*}}=\mathrm{h^{*}}/\mathrm{b^{*}}\). The  
Figure 1: Dimensionless computational domain.  
crosswind is represented by the velocity U*\(\infty\). In general, relations around h*/b* = 2.0 are typical for the study of the potential hazard to aircraft when landing or taking off between parallel runways; in that situation, large aircraft usually develop velocity of approximately U\({}_{\text{a}}\)* = 240 km/h [21].  
All quantities are nondimensionalized by adopting the length scale b*, the representative velocity \(\Gamma^{*}\)/b*, and the time scale b*/U*\(\infty\) (with crosswind) or b*\({}^{2}\)/\(\Gamma^{*}\) (with no crosswind).  
Every primary vortical structure shedding from the wingtip is discretized and represented by N Lamb discrete vortices or vortex blobs [31] randomly created there; the strength of one vortex blob is defined by \(\Gamma^{*}\)/N. Vorticity is also instantaneously generated from the ground plane. Therefore, it must be represented by M nascent vortex blobs during every time-step, \(\Delta\)t, so as to satisfy the non-slip condition in specific points (or \"pivotal points\") on the solid surface, \(\mathcal{S}_{\text{b}}\) (Figure 1). On the other hand, the method of images [32] guarantees the impermeability condition at all points on the ground plane surface. A vortex blob is completely defined in two dimensions by position vector, vortex strength, \(\Gamma_{\text{v}}\), viscous core, \(\sigma_{0}\), and vector velocity [31]. In a Lagrangian manner, every vortex blob must be individually followed during a typical numerical simulation.  
The ground plane is divided into M flat panels of the same length, namely \(\Delta\)s*, with each equipped with a centered point, the so-called \"pivotal point\". The dimensionless number \(\varepsilon/\Delta\)s (\(\varepsilon\) is the physical roughness height of the ground plane) is used to capture surface roughness effects of the ground plane in the simulation. Based on the typical roughness distribution of the ground plane, different roughness values can be chosen to be combined with different crosswind velocities, as will be discussed in Section 4.