1. Introduction
Many animals improve their locomotion energetics or increase their range by swimming or flying near the seafloor or water surface, thereby taking advantage of ground effect (Baudinette & Schmidt-Nielsen Reference Baudinette and Schmidt-Nielsen1974; Blake Reference Blake1979; Hainsworth Reference Hainsworth1988; Webb Reference Webb1993, Reference Webb2002; Park & Choi Reference Park and Choi2010). When wings/fins move steadily near a boundary, such as those of gliding animals or fixed-wing aircraft, it is well-known that a steady ground effect increases their lift and reduces their induced drag (Widnall & Barrows Reference Widnall and Barrows1970; Yang, Yang & Li Reference Yang, Yang and Li2009; Boschetti et al. Reference Boschetti, Cárdenas, Amerio and Arévalo2010; Yang, Lin & Yang Reference Yang, Lin and Yang2010; Lee & Lee Reference Lee and Lee2011; He, Qu & Agarwal Reference He, Qu and Agarwal2015; Qu et al. Reference Qu, Lu, Guo, Liu and Agarwal2015; Tremblay Dionne & Lee Reference Tremblay Dionne and Lee2018; Baddoo et al. Reference Baddoo, Kurt, Ayton and Moored2020). However, when swimmers or flyers oscillate their appendages near a boundary, like flatfish swimming near the ocean floor, they experience an unsteady ground effect.
Unsteady ground effect produces a host of hydrodynamic phenomena. For example, foils oscillating near a boundary experience a thrust or cruising speed enhancement with no penalty to their propulsive efficiency (Tandida Reference Tandida2001; Iosilevskii Reference Iosilevskii2008; Quinn, Lauder & Smits Reference Quinn, Lauder and Smits2014a; Wu et al. Reference Wu, Shu, Zhao and Yan2014; Fernández-Prats et al. Reference Fernández-Prats, Raspa, Thiria, Huera-Huarte and Godoy-Diana2015; Dai, He & Zhang Reference Dai, He and Zhang2016; Park, Kim & Sunga Reference Park, Kim and Sunga2017), and vortex wakes are deflected near boundaries due to their image vortices (Quinn et al. Reference Quinn, Moored, Dewey and Smits2014b; Kurt et al. Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019; Zhong et al. Reference Zhong, Han, Moored and Quinn2019). Moreover, stable equilibrium altitudes arise where an increase/decrease in a swimmer's altitude causes a negative/positive lift force, pulling the swimmer back to a stable altitude (Quinn et al. Reference Quinn, Moored, Dewey and Smits2014b; Mivehchi, Dahl & Licht Reference Mivehchi, Dahl and Licht2016). Kurt et al. (Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019) showed that pitching foils swimming freely in the cross-stream direction do indeed settle into a stable near-ground equilibrium altitude under no control. By showing good agreement between the equilibrium altitudes calculated through inviscid potential flow simulations and those measured in viscous water channel experiments, they determined that stable equilibrium altitudes are an inviscid phenomenon (Kurt et al. Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019). Zhong et al. (Reference Zhong, Han, Moored and Quinn2019) broadened our understanding of these equilibrium altitudes by studying hydrofoils of varying aspect ratio. They discovered that near-ground equilibrium altitudes can be stable, unstable or absent, depending upon the aspect ratio and Strouhal number.
The equilibrium altitude can determine the relative spacing between the swimmer and the boundary, which can significantly affect the hydrodynamics and propulsive performance of near-ground swimmers. To better design and control the equilibrium altitude, now, it is essential to understand the underlying physics that give rise to equilibrium altitudes. In this study, we uncover the fundamental forces that balance to generate equilibrium altitudes. Following the classic force decomposition of von Kármán & Sears (Reference von Kármán and Sears1938) and Mccune & Tavares (Reference Mccune and Tavares1993), we advance our understanding in two important ways: (i) we reveal the mechanism behind equilibrium altitudes and (ii) we leverage this understanding to develop scaling laws to predict equilibrium altitudes from only the input kinematics. Previously, Kurt et al. (Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019) proposed that stable equilibrium altitudes are generated by a balance of the added mass lift, wake-induced lift and quasi-steady lift. However, their force decomposition approach implicitly assumed that the local effective flow velocity acting on a foil was unaltered due to the presence of the ground. In this study, we no longer make this assumption and, in fact, we determine that this effect is integral to our understanding of the mechanism behind equilibrium altitudes. We will show that, upon deeper investigation, equilibrium altitudes exist due to a balance between the positive wake-induced lift and negative quasi-steady lift, with added mass lift playing precisely zero role – effectively a scallop theorem result (Purcell Reference Purcell1977) for the time-averaged lift. In this way, our new understanding of the stable equilibrium altitude mechanism supersedes that which was proposed in Kurt et al. (Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019).
2. Methods
2.1. Problem formulation
In this study, a simple model of a pitching hydrofoil is used to represent a swimmer propelling itself in an unsteady ground effect (figure 1). The foil is a 
$7\,\%$ thick teardrop-shaped cross-section (Godoy-Diana, Aider & Wesfreid Reference Godoy-Diana, Aider and Wesfreid2008; Quinn et al. Reference Quinn, Moored, Dewey and Smits2014b; Kurt et al. Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019) with a chord length of 
$c=0.095$ m. The foil undergoes sinusoidal pitching about its leading edge
\begin{equation} \theta(t)=\theta_{0}\sin(2{\rm \pi} f t), \end{equation}
where 
$\theta (t)$ is the instantaneous pitching angle, 
$\theta _{0}$ is the pitching amplitude, 
$f$ is the pitching frequency and 
$t$ is time. The pitching angle is considered positive for clockwise pitching rotations about the leading edge. The trailing-edge peak-to-peak amplitude 
$A$ is defined by
\begin{equation} A = 2c\sin\theta_{0}, \end{equation}
and its dimensionless counterpart is 
$A^* = A/c$. The reduced frequency and the Strouhal number are defined as
\begin{equation} k=\frac{fc}{U_{\infty}},\quad St = \frac{fA}{U_{\infty}}, \end{equation}
where 
$U_\infty$ is the swimming speed of the foil. The ground proximity, 
$D$, is considered as the distance from the foil's leading edge to the ground plane (figure 1) and its dimensionless form is 
$D^* = D/c$. Finally, the lift force 
$L$ is normalized by the dynamic pressure to give lift coefficient,
\begin{equation} C_{L} = \frac{L}{1/2\rho cs U_{\infty}^2}, \end{equation}
where 
$\rho$ is the fluid density while the span length 
$s$ is set to unity for two-dimensional simulations, and 
$s = {A{\kern-4pt}R} c$ for three-dimensional simulations, with 
${A{\kern-4pt}R}$ being the aspect ratio. Its time-averaged form is given as
\begin{equation} \bar{C}_{L} = \frac{\bar{L}}{1/2\rho cs U_{\infty}^2}, \end{equation}
where 
$\bar {L}$ is the time-averaged lift.
Illustration of a hydrofoil pitching near the ground. An inertial frame of reference fixed to the undisturbed fluid is denoted by 
$(X,Z)$ while a body-fixed frame of reference is denoted by 
$(x,z)$. An image system enforces a no-flux boundary condition at the ground plane (
$Z = 0$).
2.2. Numerical methods
The simulations model a potential flow that is assumed to be irrotational, incompressible and inviscid. This flow is governed by Laplace's equation, 
$\nabla ^{2}\varPhi ^{*}=0$, where 
$\varPhi ^{*}$ is the perturbation velocity potential in an inertial frame fixed to the undisturbed fluid. A no-flux boundary condition, 
$\boldsymbol {\nabla } \varPhi ^*\boldsymbol{\cdot }\boldsymbol {n}=0$, where 
$\boldsymbol {n}$ is the surface normal vector, must be satisfied on the body surface and the ground plane at 
$Z=0$, and a far-field boundary condition such that flow perturbations decay far away from the body and its wake surfaces must be satisfied.
Following previous studies (Katz Reference Katz1985; Quinn et al. Reference Quinn, Moored, Dewey and Smits2014b; Moored Reference Moored2018; Kurt et al. Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019; Zhong et al. Reference Zhong, Han, Moored and Quinn2019; Mivehchi et al. Reference Mivehchi, Zhong, Kurt, Quinn and Moored2021), this potential flow problem is solved numerically using an unsteady boundary element method with the method of images used to enforce the ground plane boundary condition. In this method, the foil, its wake surfaces and their images are discretized into a finite number of boundary elements. The general solution to the potential flow problem is then reduced to finding a distribution of constant-strength doublet and source elements that satisfy the body and ground plane boundary conditions. The boundary elements also implicitly satisfy the far-field boundary condition. The body boundary condition is imposed using a Dirichlet condition at collocation points centred on each body element, but pushed within the body along the element's surface normal. Also, an explicit Kutta condition is applied where the vorticity at the trailing edge is zero by tuning the strength of a wake boundary element at the trailing edge. At each time step, the trailing-edge element is shed with a strength that satisfies Kelvin's circulation theorem and it remains constant thereafter. A wake roll-up algorithm is implemented to advect the wake elements with the local velocity calculated by the desingularized Biot–Savart law (Krasny Reference Krasny1986; Zhu et al. Reference Zhu, Wolfgang, Yue and Triantafyllou2002; Moored Reference Moored2018; Akoz, Mivehchi & Moored Reference Akoz, Mivehchi and Moored2021). Therefore, a matrix representation of the boundary condition is formed and solved for the body doublet strengths at each time step. The perturbation velocity over the surface of the body is then determined by a local differentiation of the perturbation potential.
Following the work of Moored (Reference Moored2018), the pressure exerted on the body is calculated by the unsteady Bernoulli equation as
\begin{equation} P(x,z,t)=-\rho \frac{\partial \varPhi^*}{\partial t} + \rho (\boldsymbol{u_{rel}} + \boldsymbol{U_{\infty}}) \boldsymbol{\cdot } \boldsymbol{\nabla} \varPhi^* -\rho \frac{(\nabla \varPhi^*)^2}{2}, \end{equation}
where 
$\boldsymbol {u_{rel}}$ is the local surface velocity of the pitching foil relative to the body-fixed frame of reference. The inviscid instantaneous resultant force is then simply an integration of the instantaneous pressure field acting over the body surface. Finally, the resultant force is decomposed into its lift and drag directions. For more details see Moored (Reference Moored2018).
2.3. Force decomposition approach
Following the classic force decomposition of von Kármán & Sears (Reference von Kármán and Sears1938) and Mccune & Tavares (Reference Mccune and Tavares1993), the lift force can be decomposed as the sum of the added mass lift, quasi-steady lift and wake-induced lift, which are denoted by 
$({\cdot })_{a}$, 
$({\cdot })_{q}$ and 
$({\cdot })_{w}$, respectively. We apply this decomposition approach directly to the pressure field defined by the unsteady Bernoulli equation.
To calculate the added mass pressure, 
$P_a$, potential flow simulations are employed where there are no wake elements shed and no Kutta condition enforced. Under these conditions, the first term in (2.6) is known as the added mass pressure. To calculate the quasi-steady pressure, 
$P_q$, potential flow simulations are employed where there are no wake elements shed, but there is a Kutta condition enforced (modelled in the boundary element method with a trailing-edge panel that extends far from the foil, effectively to infinity). Then the combined second and third terms in (2.6) provide the quasi-steady pressure. The combined added mass and quasi-steady pressures are
\begin{align} P_a(x,z,t) + P_q(x,z,t)=\left[\underbrace{-\rho \frac{\partial \varPhi^*}{\partial t}}_{{added\ mass}} + \underbrace{\rho (\boldsymbol{u_{rel}} + \boldsymbol{U_{\infty}}) \boldsymbol{\cdot } \boldsymbol{\nabla} \varPhi^* -\rho \frac{(\nabla \varPhi^*)^2}{2}}_{quasi\text-steady}\right]_{{NO\ WAKE}}. \end{align}
The wake-induced pressure, 
$P_w$, is calculated as the difference between the total pressure, 
$P$ (calculated from the full simulations described in § 2 by (2.6)), and the combined added mass and quasi-steady pressures
\begin{equation} P_w(x,z,t)=P(x,z,t) - \left[P_a(x,z,t) + P_q(x,z,t)\right]. \end{equation}Finally, each force component is acquired by integrating its corresponding pressure field over the body. Then the added mass, quasi-steady and wake-induced lift coefficients can be defined and they sum to the total lift coefficient
\begin{equation} C_{L} = C_{L_{a}} + C_{L_{q}} + C_{L_{w}}. \end{equation}2.4. Force decomposition validation
Before applying the force decomposition approach to simulations of a pitching foil in ground effect, we first validated it against Theodorsen's solution, which provides a test case where the added mass, quasi-steady and wake-induced terms are readily identified (Bisplinghoff, Ashley & Halfman Reference Bisplinghoff, Ashley and Halfman2013). In Theodorsen (Reference Theodorsen1935), the lift force generated by a foil pitching about its leading edge is
\begin{align} L &= \underbrace{\rho {\rm \pi}\frac{c^2}{4}\left(U_{\infty} \dot{\theta} + \frac{c}{2} \ddot{\theta}\right)}_{{added\ mass}} + \underbrace{\rho {\rm \pi}U_{\infty} c \left(U_{\infty}\theta + \frac{3c}{4}\dot{\theta}\right)}_{quasi\text-steady} \nonumber\\ &\quad +\underbrace{\left[C(k)-1\right]\rho {\rm \pi}U_{\infty} c \left(U_{\infty}\theta + \frac{3c}{4}\dot{\theta}\right)}_{wake\text-induced}, \end{align}
where 
$C(k)$ is the lift deficiency function. Full boundary element simulations and their decomposition counterparts were performed as described in §§ 2.2 and 2.3 on a thin airfoil (
$1\,\%$ thick teardrop shape), undergoing small pitching amplitudes of motion about its leading edge (
$A^* =0.002$), with a frozen wake (no wake roll-up algorithm applied), and across a wide reduced frequency range of 
$k = 0.1, 1$ and 
$10$. Figure 2 shows that all three lift components precisely match the Theodorsen solution across a wide range of 
$k$. Even though Theodorsen's solution makes several assumptions, it is worth noting that the force decomposition approach itself is generic and has no underlying assumptions.
Validation of the force decomposition approach against Theodorsen's solution across a wide reduced frequency range of 
$0.1 \leq k \leq 10$. The panel rows show: (a) added mass lift, (b) wake-induced lift and (c) quasi-steady lift.
2.5. Experimental methods
We conducted four groups of experiments to facilitate the force decomposition analysis, including a two-dimensional (2-D) static ground effect test, a 2-D unsteady high-frequency test (for wake deflection analysis) and a 2-D and a 3-D unsteady ground effect test. All hydrofoils had a teardrop cross-section with a chord length of 
$0.095$ m to ensure consistency with simulations, and were 3D-printed with solid ABS (Acrylonitrile Butadiene Styrene). A stainless steel driveshaft with a diameter of 
$6.35$ mm was used to pitch the hydrofoil at its leading edge, with the pitch angle 
$\theta$ being prescribed as 
$\theta _0 \sin (2{\rm \pi} f t)$. The driveshaft was driven by a high-torque digital servo motor (Dynamixel MX64), while an absolute encoder (US Digital A2K 4096 CPR) was utilized to verify the angle of the actual motion. The forces and moments generated were measured using a six-axis force-torque sensor (ATI-Mini 40: SI-40-2). The lateral forces were time averaged over 30 seconds or 20 oscillation cycles to determine the mean lift force for static and unsteady tests. The incoming flow speed was set to 143 mm s
$^{-1}$ using an ultrasonic flowmeter (Dynasonics Series TFXB), resulting in a chord-based Reynolds number of 13 500. For ground effect tests, an additional splitter plate was used instead of a tunnel sidewall to control the boundary layer thickness (figure 3). A boundary layer thickness of 
$\delta _{99\,\%} \approx 7.5$ mm (
$\delta _{99\,\%}/c = 0.08$) was measured using particle image velocimetry at the position aligned with the leading edge of the hydrofoil.
Schematic of experiment set-up equipped with removable splitter plate.
For all 2-D tests, we used a hydrofoil of aspect ratio 3 in a nominally 2-D flow condition, which was achieved by installing a horizontal splitter plate and a surface plate near the tips of the hydrofoil. The gap between the hydrofoil tips and the surface/splitter plate was less than 5 mm and surface waves were also minimized by the presence of the surface plate.
To validate the 3-D force decomposition, we repurposed the 3-D unsteady ground effect data that were previously published in Zhong et al. (Reference Zhong, Han, Moored and Quinn2019). We tested hydrofoils with four aspect ratios: 1, 1.5, 2 and 2.5, and assessed their performance over a range of ground proximities. Each trial was repeated three times.
3. Force decomposition results
A pitching foil in ground effect and constrained to fixed ground distances was simulated across a wide range of dimensionless ground proximities, reduced frequencies and Strouhal numbers, and force decomposition was applied as described in § 2.3. The simulation variables are detailed in table 1.
Simulation variables investigated in the current study.
3.1. Mechanism behind the generation of stable and unstable equilibria
The pitching foil examined in this study is effectively out of ground effect by 
$D^* = 1.5$ for cases with 
$St<0.45$, and by 
$D^* = 2$ for cases with higher 
$St$. In fact, a typical reverse von Kármán vortex street, indicative of thrust production in an infinite domain, is observed when 
$D^* = 1.5$, 
$k = 1$ and 
$St = 0.3$ (figure 4a). As the foil approaches the ground to a distance of 
$D^* = 0.4$ (
$k = 1$ and 
$St = 0.3$), the wake vortices deflect away from the ground (figure 4b), an effect that has been attributed to the image vortex system (Quinn et al. Reference Quinn, Moored, Dewey and Smits2014b).
Both the stable and unstable equilibrium altitudes are generated by a balance between the positive wake-induced lift and negative quasi-steady lift. (a) Wake plots for an out-of-ground effect foil. (b) Wake plots for a near-ground foil. (c) Total lift for varying 
$k$ with 
$St =0.3$. (e) Total lift for varying 
$St$ with 
$k = 1$. (d) Zoom-in plot highlighting the unstable equilibrium altitude. (f,i) Added mass lift. (g,j) Wake-induced lift. (h,k) Quasi-steady lift. Note that the stable equilibrium altitude is marked by a solid line, the stable equilibrium altitude for a foil unconstrained in the cross-stream direction (mass of the foil is chosen to be 
$m=2.68\rho s c^2$) is marked by the dashed line and the unstable equilibrium altitude is marked by the dotted line.
The total lift for all of the simulations is presented for varying reduced frequency and ground distance (
$St = 0.3$; figure 4c), and for varying Strouhal number and ground distance (
$k = 1$; figure 4d,e). A common trend emerges in all of the in-ground effect data except when 
$St > 0.475$: the lift is positive close to the ground (pushing a foil away from the ground) and negative far from the ground (pulling a foil towards the ground), crossing zero (stable equilibrium altitude; solid lines in figure 4) between these two regimes. As observed previously, the stable equilibrium altitude increases with increasing 
$St$ and decreasing 
$k$ (Kurt et al. Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019), and an unstable equilibrium altitude emerges when 
$St > 0.475$ (dotted line in figure 4d), which increases slightly with increasing 
$St$ in a narrow near-ground region (Zhong et al. Reference Zhong, Han, Moored and Quinn2019). Altitudes below/above the unstable equilibria produce a negative/positive lift, demonstrating that it is indeed unstable. For the first time, this unstable equilibrium altitude is observed in simulation, which corroborates previous experimental findings (Zhong et al. Reference Zhong, Han, Moored and Quinn2019).
The dashed lines in figures 4(c) and 4(e) show the unconstrained simulation results from Kurt et al. (Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019). As discovered in Cochran-Carney (Reference Cochran-Carney2018), the observed difference in the equilibrium altitudes of the constrained and unconstrained cases is ultimately caused by a difference in the unconstrained simulation's dimensionless mass (
$m^* \equiv m/(\rho s c^2) = 2.68$) compared with that of the constrained simulations (
$m^* = \infty$). The finite mass foils of the unconstrained simulations have a finite heaving recoil while the effectively infinite mass foils of the constrained simulations have zero heaving recoil, thereby altering the force production between the two cases and, consequently, their equilibrium altitudes. Despite these differences, the constrained equilibria show all of the same trends as the unconstrained equilibria with only a small shift in their 
$D^*$ location. Moreover, this study will focus on constrained simulations by design since these allow us to investigate the forces outside of equilibrium conditions, which is essential to fully understand the stability properties of the equilibria.
The decomposition of the lift coefficient into its added mass, wake-induced and quasi-steady components (figure 4f–k) shows clear trends in the components for all ground proximities and kinematics. The added mass lift is always precisely zero and therefore does not play a role in generating equilibrium altitudes (figure 4f,i). For in-ground effect data, the wake-induced lift is always positive (figure 4g,j) while the quasi-steady lift is always negative (figure 4h,k), revealing that both unstable and stable equilibria are generated by a balance between positive wake-induced lift and negative quasi-steady lift. The magnitudes of the wake-induced and quasi-steady lift increase with increasing 
$St$ and decreasing 
$k$.
The data contours in figure 4 reveal the basic mechanism behind the existence of equilibria, but they do not provide a clear picture as to why there can be both a stable and an unstable equilibria for a given 
$St$ and 
$k$. To better understand this feature of the data, the lowest and highest 
$St$ cases with 
$k = 1$, 
$St = 0.25$ and 
$St = 0.55$, respectively, are presented in figure 5. For low 
$St$ (figure 5a), the total lift only crosses zero at one 
$D^*$, producing a stable equilibrium. This equilibrium altitude is stable since, for altitudes above equilibrium, the magnitude of the quasi-steady lift is greater than the wake-induced lift, leading to a net negative total lift force, and for altitudes below equilibrium the magnitude of the wake-induced lift is greater than the quasi-steady lift, leading to net positive total lift. For high 
$St$ (figure 5b), the total lift is seen to cross through two equilibrium altitudes. Now, it becomes clear that the wake-induced lift magnitude curve crosses through the quasi-steady lift magnitude curve at two ground distances because the wake-induced lift has a nearly linear increase with decreasing ground distance while the quasi-steady lift has an exponential growth for 
$D^*<1.1$. In fact, the slopes of the lift components at the equilibrium points determine their stability: when 
$|{\rm d}\bar {C}_{Lq}/{\rm d} D^*| < |{\rm d}\bar {C}_{Lw}/{\rm d} D^*|$ the equilibrium altitude is stable and when 
$|{\rm d}\bar {C}_{Lq}/{\rm d} D^*| > |{\rm d}\bar {C}_{Lw}/{\rm d} D^*|$ the equilibrium altitude is unstable.
Force decomposition for 
$k = 1$: (a) 
$St = 0.25$ and (b) 
$St=0.55$.
Although this analysis has revealed the basic mechanism behind the generation of stable and unstable equilibrium altitudes, it also raises several questions. How does added mass lift integrate precisely to zero when we expect peak negative added mass lift to be larger in magnitude than the peak positive added mass lift since the foil will get closer and further from the ground throughout an oscillation cycle? Why is the wake-induced lift positive for a wake deflected away from a ground plane when, based on momentum flux arguments, we would expect the opposite? Why is the quasi-steady lift negative in nature? To address these questions, the following sections will examine each force component more deeply.
3.2. Added mass lift
As shown in the previous section, the time-averaged added mass lift is precisely zero for all of the simulation cases. This is expected for solutions to the linear potential flow equations with harmonic boundary conditions – a ‘scallop theorem’ result (Purcell Reference Purcell1977). However, for the in-ground effect data it is well-known that, as the ground distance decreases, there is an increase in the added mass of a foil (Brennen Reference Brennen1982; Mivehchi et al. Reference Mivehchi, Zhong, Kurt, Quinn and Moored2021). This effect should amplify both the positive and negative peaks in the added mass lift, but should more greatly amplify the negative peak, which occurs at the bottom of the down stroke when the foil is closest to the ground. Considering only these peak forces, one would imagine that the time-averaged added mass lift would be negative, yet we know that it is zero from the data in § 3.1, so how does the added mass lift integrate precisely to zero and satisfy a scallop theorem?
Figure 6 presents the time-varying added mass lift at four ground distances when 
$k = 1$ and 
$St = 0.3$. Indeed, as the ground distance decreases, both the positive and negative peaks grow in magnitude, and it is true that the negative peak is amplified more than the positive peak, creating a net negative lift in the peak forces. However, it is now clear that the negative peak becomes sharper while the positive peak becomes broader. When integrating the lift over a full cycle, this effect counteracts the net imbalance in the positive and negative peak forces, leading to a precisely zero time-averaged force.
The added mass lift as a function of time for varying ground distance with 
$k = 1$ and 
$St = 0.3$. Out-of-ground effect data are represented by 
$D^* = \infty$.
Although our simulations were for potential flow, the phenomenon of the time-averaged added mass lift being precisely zero is also valid in viscous flow. According to Graham, Pitt Ford & Babinsky (Reference Graham, Pitt Ford and Babinsky2017), in viscous flow, the added mass lift of a foil is the time derivative of the quasi-steady vortex impulse in the cross-stream direction, which is
\begin{equation} L_{a} =-\rho\frac{\text{d}I_{z}}{\text{d}t},\quad \text{and}\quad I_{z} =-\int x \gamma_{0} \,\text{d}S, \end{equation}
where 
$x$ is the coordinate of the body-fixed frame of reference, 
$\gamma _{0}$ is the bound quasi-steady vortex sheet, 
$S$ is the body surface and 
$I_{z}$ is the vortex impulse. These equations suggest that, since the quasi-steady bound circulation is a periodic signal for a foil periodically pitching near the ground, then the vortex impulse is also a periodic signal. When the added mass lift is then time averaged over a period of motion, it is necessarily zero since it is the time derivative of the vortex impulse.
3.3. Wake-induced lift
In § 3.1, the wake-induced lift was found to be positive for in-ground effect foils regardless of their ground distance and kinematics. This is surprising, since the wake vortices deflect away from the ground and in the time average they also produce a momentum jet deflected away from the ground (Quinn et al. Reference Quinn, Moored, Dewey and Smits2014b). If one imagines a control volume around a foil and considers only the momentum flux term, the deflection of the jet away from the ground should have an associated negative lift force, that is, a force pushing the swimmer toward the ground (Kurt et al. Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019), and yet the data from § 3.1 show otherwise.
To understand this phenomenon more deeply we employ a control volume analysis. To begin, we will consider a simple out-of-ground effect case where a deflected wake is also generated very much akin to the deflected wake of in-ground effect foils. This case is useful since it will help unambiguously show the connection among the wake-induced lift, the pressure field and the momentum flux, and it presents an opportunity to experimentally verify the trends in the wake-induced lift. With this basic knowledge in mind, we will then proceed to the in-ground effect case.
When a 2-D pitching foil oscillates at sufficiently high Strouhal number, it is well documented (Marais et al. Reference Marais, Thiria, Wesfreid and Godoy-Diana2012; Das, Shukla & Govardhan Reference Das, Shukla and Govardhan2016) that there is a wake instability that causes the wake vortices to form vortex dipoles and, in turn, deflect the wake. Figure 7(a) shows an example of this, where a pitching foil is operated at 
$St = 0.55$ and 
$A^* = 0.3$. Since this pitching foil is out of ground effect and has a symmetric motion, then the time-averaged added mass and quasi-steady lift must both necessarily be zero. However, due to the wake deflection the time-averaged wake-induced lift will be non-zero and, in fact, it will be the only component thus making it equivalent to the total lift in this specific case.
Time-averaged lift control volume analysis. The control volume is represented by a dashed box, where 
$C_{z}$ and 
$C_{x}$ are the vertical and horizontal dimensions of the volume, respectively. Here, 
$C_{z}$ is fixed at 
$2c$ and 
$C_{x}$ is fixed at 
$4c$ in (c,d), but varied in (e,f). (a) Upward deflected wake generated by an out-of-ground effect foil pitching at 
$St = 0.55$ and 
$A^*=0.3$. (b) Upward deflected wake generated by an in-ground effect foil pitching at 
$St = 0.3$ and 
$A^*=0.3$. (c) Control volume analysis for the out-of-ground effect case. (d) Control volume analysis for the in-ground effect case. (e) Size effect of the control volume for the in-ground effect case. (f) Size effect of the control volume for the out-of-ground effect case. The time-averaged wake-induced lift coefficient from simulations and experiments are denoted by the blue lines and blue circles, respectively. Contributions to simulations’ wake-induced lift from the pressure difference term and the momentum flux term in the control volume analysis are shown as green and yellow lines, respectively.
Following the work of Cleaver, Wang & Gursul (Reference Cleaver, Wang and Gursul2012), a time-averaged control volume analysis can be applied to the control volume (dashed box) defined in figure 7(a)
\begin{equation} \bar{L}_{w}=\underbrace{\overline{\int_{UL} \Delta P_{UL}\,{\rm d} A}}_{{Pressure\ difference}} \underbrace{-\overline{\int_{CS} \rho v(\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{n})\,{\rm d} A}}_{{Momentum\ flux}}, \end{equation}
where the total time-averaged lift is precisely equal to the time-averaged wake-induced lift, 
$\bar {L}_{w}$, the pressure difference between the upper and lower (
$UL$) control surfaces is 
$\Delta P_{UL}$, the velocity field is 
$\boldsymbol {V}$, the vector normal to the control surfaces (
$CS$) is 
$\boldsymbol {n}$ and the horizontal and vertical components of the velocity field are 
$u$ and 
$v$, respectively. In this sense, the wake-induced lift consists of a momentum flux term and a pressure difference term. The left bottom corner of the control volume is 
$c$ from the foil's leading edge vertically and 
$0.25c$ horizontally. The vertical and horizontal dimensions of the volume, 
$C_{z}$ and 
$C_{x}$, are fixed at 
$2c$ and 
$4c$, respectively.
Just like the in-ground effect data from § 3.1, both simulations and experiments show that the wake-induced lift is in fact positive for an upward deflected wake (figure 7c) which, in the case of this out-of-ground effect foil, is consistent with previous findings (Emblemsvåg, Suzukiy & Candlerz Reference Emblemsvåg, Suzukiy and Candlerz2002; Liang et al. Reference Liang, Ou, Premasuthan, Jameson and Wang2011; Cleaver et al. Reference Cleaver, Wang and Gursul2012; Yu, Hu & Wang Reference Yu, Hu and Wang2012). Using the simulation data, figure 7(c) further shows that the wake-induced lift from the net momentum flux term is indeed negative, as expected, while the contribution from the pressure difference term is positive and, surprisingly, is larger in magnitude than the momentum flux term, thereby making the wake-induced lift positive for an upward deflected wake. In light of this analysis, momentum flux considerations alone do not accurately discern the direction of forces generated by an asymmetric wake; the pressure field must also be considered.
Moreover, the same control volume analysis was applied to an in-ground effect pitching foil in figure 7(b), which was calculated by (3.3) and (3.4). The left bottom corner of the control volume is 
$D$ from the foil's leading edge vertically and 
$0.25c$ horizontally. The vertical and horizontal dimensions of the volume, 
$C_{z}$ and 
$C_{x}$, are fixed at 
$2c$ and 
$4c$, respectively. Since the time-averaged added mass lift is zero (§ 3.2), the time-averaged wake-induced lift from the pressure difference and momentum flux terms are acquired by calculating those two contributions for the total and quasi-steady pressure and flow fields and taking their differences as
\begin{equation} \bar{L}_{w,{pressure}}= \left.\overline{\int_{UL} \Delta P_{UL}\,{\rm d} A}\right|_{{total}} - \left.\overline{\int_{UL} \Delta P_{UL}\,{\rm d} A}\right|_{quasi\text-steady}, \end{equation}and
\begin{equation} \bar{L}_{w,{momentum}} = \left.-\overline{\int_{CS} \rho v(\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{n})\,{\rm d} A}\right|_{{total}} + \left.\overline{\int_{CS} \rho v(\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{n}),{\rm d} A}\right|_{quasi\text-steady}. \end{equation}Figure 7(d) shows that, in ground effect, the momentum flux term is negative, yet the wake-induced lift acts in the positive direction since the pressure difference term again outweighs the momentum flux term, as in the out-of-ground effect case. It is worth noting that compared with the out-of-ground effect case the pressure difference term now has a more significant effect on the wake-induced lift.
To show that the conclusion of the pressure difference term determining the direction of the wake-induced lift is not affected by the size of the control volume, 
$C_{x}$ of the control volume is varied from 
$3.5c$ to 
$4.5c$ for both out-of-ground and in-ground cases numerically (figure 7e,f). It is discovered that the magnitudes of the two terms stay nearly constant, and still, the magnitude of the positive pressure difference term is larger than that of the negative momentum flux term. Our simulations and the experiments in Cleaver et al. (Reference Cleaver, Wang and Gursul2012) suggest that, as long as the momentum flux crosses the control volume boundary, the conclusion that the pressure difference term determines the direction of the wake-induced lift still holds.
3.4. Quasi-steady lift
To provide deeper insight as to why there is a net negative quasi-steady force in the time average, we further break down the quasi-steady lift into its constituent components by leveraging the Kutta–Joukowski theorem (von Kármán & Sears Reference von Kármán and Sears1938), 
$L_{q} = \rho s U \varGamma _{0}$, where 
$L_q$ is the quasi-steady lift, 
$s$ is the span length, 
$U$ is the local effective flow velocity acting on a foil and 
$\varGamma _{0}$ is the quasi-steady bound circulation of a foil. The local effective velocity is calculated from simulations by summing the free-stream velocity, 
$U_\infty$, with the induced velocity 
$U_{i}$ at the leading edge from the image body. The quasi-steady bound circulation is simply the negative of the strength of the trailing-edge panel used to enforce the Kutta condition in the boundary element simulations, that is, 
$\varGamma _0 = -\mu _{TE}$ (Katz & Plotkin Reference Katz and Plotkin2005). The Kutta–Joukowski theorem can be non-dimensionalized and written as
\begin{equation} C_{L_q} = 2 U^* \varGamma_0^*,\quad U^* = \frac{U}{U_\infty},\quad \varGamma_0^* = \frac{\varGamma_0}{cU_\infty}. \end{equation}
When a foil is out of ground effect then 
$U^* = 1$, 
$\varGamma _{0}^* = \varGamma _{0,\infty }^*$ and 
$C_{L_q} = C_{L_{q,\infty }}$, where 
$\varGamma _{0,\infty }^*$ and 
$C_{L_{q,\infty }}$ are the time-varying quasi-steady bound circulation and lift coefficient for a foil in an infinite domain. The quasi-steady bound circulation has been found to be amplified as the foil gets closer to the ground (Zong, Liang & Zhou Reference Zong, Liang and Zhou2012; Baddoo et al. Reference Baddoo, Kurt, Ayton and Moored2020). This effect comes from the fact that, as the foil approaches the ground, the bound vortex of its image foil induces a transverse velocity at the trailing edge of the real foil. In order for the Kutta condition to be preserved, an additional bound circulation of the real foil must arise to cancel the additional transverse velocity induced by the image foil. Without loss of generality, the quasi-steady bound circulation for in-ground effect foils can then be modelled as the multiplication of a circulation amplification factor 
$\beta$ and the infinite domain quasi-steady circulation, 
$\varGamma _{0}^* = \beta \varGamma _{0,\infty }^*$. We expect that both the circulation and the effective flow velocity will be amplified or reduced for foils in ground effect as opposed to out-of-ground effect foils.
To understand the connections among the circulation amplification factor, the amplified/reduced effective flow velocity and the quasi-steady lift, as well as to verify these connections with experimental data, we begin with a simple case of a static foil in and out of ground effect at various pitch angles (figure 8). For the static foil, the quasi-steady lift is precisely the total lift, which is readily measured in experiments. Figure 8(a) presents the quasi-steady lift for 
$D^* = 0.3$ and 
$D^* = \infty$ from the simulations (solid lines), which is directly determined by integrating the pressure acting on the foil calculated from the unsteady Bernoulli equation (§ 2.2). The simulations are also used to calculate the quasi-steady lift by using the Kutta–Joukowski theorem (circle markers) as outlined above, and both approaches are observed to produce exactly the same result. Also, the in-ground effect quasi-steady lift is verified by having good agreement with experimental data (square markers).
Asymmetry in the quasi-steady lift for a steady foil in ground effect is caused by the velocity induced on a foil by its image body's quasi-steady bound circulation. (a) The quasi-steady lift coefficient. Note that the lift coefficient calculated in the simulations (using the unsteady Bernoulli equation) is represented with solid lines, the lift coefficient calculated from the simulations using the Kutta–Joukowski theorem is represented by the circle markers and the experimental data are represented by the square markers. (b) The dimensionless quasi-steady circulation. (c) The dimensionless effective flow velocity. (d) Schematic of the components of the effective flow velocity, U, that is the summation of 
$U_{i}$ and 
$U_{\infty }$.
For the out-of-ground effect foil, the quasi-steady lift shows the classic linear relationship with the pitch angle for 
$-8^\circ \leq \theta \leq 8^\circ$. However, when in ground effect the quasi-steady lift shows a nonlinear relationship with the pitch angle, where both the negative and positive lift regimes are amplified compared with the out-of-ground effect case, yet the negative lift regime is amplified more than the positive lift regime. This produces an asymmetric lift response about 
$\theta = 0$. If one were to pitch the foil through an 
$8^\circ$ amplitude sinusoidal pitching motion at an infinitesimal Strouhal number, this asymmetric lift response would lead to a net negative quasi-steady lift in the time average (
$\bar {C}_{L_{q}}=-0.07$ for simulations and 
$\bar {C}_{L_{q}}=-0.09$ for experiments), just as is observed in the quasi-steady data of finite Strouhal numbers. In fact, throughout this section, we will show that precisely the same mechanism is at play in this static case as is in the quasi-steady cases of finite Strouhal number.
Figure 8(b) presents the dimensionless quasi-steady circulation. The circulation amplification is observed to be symmetric about 
$\theta = 0$ and appears to just be a constant factor giving rise to a slope change in the 
$\varGamma _0^*$ – 
$\theta$ line. This shows, at least for the static case, that 
$\beta$ is only a function of 
$D^*$, not 
$\theta$. Moreover, the circulation amplification is responsible for amplifying the quasi-steady lift forces in ground effect, but is not responsible for asymmetry in the lift curve. Figure 8(c) shows that the effective flow velocity increases above the free-stream speed when the quasi-steady circulation is negative (counterclockwise), while it decreases below free-stream speed when the quasi-steady circulation is positive (clockwise). This occurs due to the velocity induced by the quasi-steady bound vortex of the image body (figure 8d). By consequence, increased effective flow speed acts to amplify the negative lift of negative pitch angles (negative bound circulation) thereby acting in concert with the circulation amplification effect, and decreased effective flow speed acts to reduce the positive lift of positive pitch angles (positive bound circulation), thereby counteracting the circulation amplification effect. Therefore, it is the change in the effective flow speed that leads to the asymmetric lift response. This mechanism will be implicated further as the reason for the net negative quasi-steady lift for oscillating foils from zero to finite Strouhal numbers.
Figure 9 presents the quasi-steady lift coefficient, dimensionless quasi-steady circulation and dimensionless effective flow velocity for a foil at a range of ground distances with 
$St = 0.3$. In figure 9(a), first, the quasi-steady lift calculated by the unsteady Bernoulli equation and by the Kutta–Joukowski theorem still show good agreement. Figure 9(b) shows that the negative and positive quasi-steady circulations are symmetrically amplified with decreasing ground distance, precisely as in the static foil case. Moreover, figure 9(c) shows that the effective flow velocity increases above and decreases below the free-stream velocity when the bound circulation is negative and positive, respectively, also mirroring the static foil case. Consequently, the change in effective flow velocity leads to an enhancement in the negative lift and a reduction in the positive lift thereby generating a net negative quasi-steady lift in the time average over one pitching cycle. Figure 9(d) illustrates precisely the same effective flow altering mechanism of the image body's quasi-steady circulation as is in the static foil case.
The velocity induced by the ground image body creates asymmetry in the quasi-steady lift at finite Strouhal numbers just like the steady foil case. (a) The quasi-steady lift coefficient. (b) The dimensionless quasi-steady circulation. (c) The dimensionless effective flow velocity. (d) Schematic showing the effective flow velocity, U, is the summation of 
$U_{i}$ and 
$U_{\infty }$.
3.5. Three-dimensional decomposition
In figure 10(a), both 3-D simulations and experiments observe that the stable equilibrium altitude disappears when the aspect ratio decreases to 
${A{\kern-4pt}R} = 1$, which is reproduced from Zhong et al. (Reference Zhong, Han, Moored and Quinn2019). Figure 10(b), presents the lift decomposition of the 3-D simulations into the added mass, quasi-steady and wake-induced components. Exactly as in the 2-D simulations, a stable equilibrium altitude arises for foils of finite 
${A{\kern-4pt}R}$ by a balance between the positive wake-induced lift and negative quasi-steady lift. The equilibrium altitude disappears for low aspect ratio (
${A{\kern-4pt}R} < 1.5$) hydrofoils, which instead experience a net negative time-averaged lift force. The decomposition reveals that this occurs since the magnitude of the wake-induced lift reduces more with decreasing 
${A{\kern-4pt}R}$ than the reduction in the magnitude of quasi-steady lift, thereby having the negative quasi-steady lift outweigh the positive wake-induced lift at all ground distances.
The magnitude of the wake-induced lift of 3-D foils degrades more with decreasing 
${A{\kern-4pt}R}$ than their quasi-steady lift leading to the disappearance of an equilibrium altitude for 
${A{\kern-4pt}R} < 1.5$. Simulation data of pitching 3-D hydrofoils of varying aspect ratio with 
$St= 0.25$ and 
$\theta _{0} = 11^{\circ }$. (a) Total lift as calculated from the simulations. Note that the circle markers represent the stable equilibrium altitudes measured by experiments. (b) Decomposed lift components from the simulations.
3.6. Scaling the equilibrium altitude
The analysis of the added mass, quasi-steady and wake-induced forces presented above can now be leveraged to develop scaling laws that capture the physics of ground effect and are capable of predicting the equilibrium altitude of a pitching foil. Since the added mass forces do not play a role in the generation of an equilibrium altitude then scaling laws for only the quasi-steady and wake-induced lift are sought, two flow models were used to develop the scaling laws for the quasi-steady lift and wake-induced lift, respectively, which enables us to pinpoint any sources of the modelling error in the scaling laws.
3.6.1. Quasi-steady lift
As presented in § 3.4, the quasi-steady lift coefficient is
$$\begin{gather} C_{L_{q}} \propto U^* \beta \varGamma_{0,\infty}^*, \end{gather}$$
$$\begin{gather}\varGamma_{0,\infty}^* = {\rm \pi}\theta +\frac{3{\rm \pi} c \dot{\theta}}{4U_{\infty}}, \end{gather}$$
where 
$\beta$ is the quasi-steady circulation amplification factor, 
$U^*$ is the normalized effective flow velocity at the foil's leading edge and 
$\varGamma ^*_{0,\infty }$ is simply the formula for a flat plate pitching about its leading edge from Theodorsen (Reference Theodorsen1935) and Mccune & Tavares (Reference Mccune and Tavares1993). To determine a scaling law for the quasi-steady lift, consider a simplified flow model where the quasi-steady bound vortex is located along the chord (figure 9d). Then a scaling relation for the circulation amplification factor and the effective flow velocity can be determined. The circulation amplification factor is calculated by finding the circulation amplification necessary to enforce the Kutta condition at the foil's trailing edge due to the influence of the image quasi-steady bound vortex. Assuming that the distance between the bound vortex and foil's leading edge is 
$bc$ with 
$0\leq b \leq 1$, then the circulation amplification factor is simply calculated as
\begin{equation} \beta = 1 + \frac{(1-b)^2}{4{D^*}^2} \quad \text{where}\ b = 0.55. \end{equation}
To determine the placement of the bound vortex along the chord, circulation amplification data at 
$D^*=0.5$ from an exact solution for a steady foil in ground effect (figure 4(d) in Baddoo et al. Reference Baddoo, Kurt, Ayton and Moored2020) were used to determine that 
$b=0.55$. The effective flow velocity is calculated by finding the normalized induced velocity at the leading edge, 
$U^*_{i} =U_{i}/U_\infty$, from the image foil's quasi-steady bound vortex
\begin{equation} U^* = 1 + U^*_{i} \quad \text{where}\ U^*_{i} = \frac{-\varGamma_{0,\infty}^*\beta}{\rm \pi}\left(\frac{D^*}{b^2 + 4{D^*}^2}\right). \end{equation}After organization, the time-varying quasi-steady lift becomes
\begin{equation} C_{L_{q}} \propto \beta\varGamma_{0,\infty}^* - \frac{\beta^2{\varGamma^*}^2_{0,\infty}D^*}{{\rm \pi}(b^2+4{D^*}^2)}. \end{equation}
According to thin airfoil theory, the coefficient of the first term is 
$2$, which gives
\begin{equation} C_{L_{q}} = 2 \beta\varGamma_{0,\infty}^*-c_{1} \frac{\beta^2{\varGamma^*}^2_{0,\infty}D^*}{{\rm \pi}(b^2+4{D^*}^2)},\end{equation}
where 
$c_{1}$ needs to be determined.
As an important side note, this formula can be used to determine a lift scaling relation for steady ground effect by neglecting the rotational lift or virtual camber term (
$\dot {\theta }$ term) in 
$\varGamma ^*_{0,\infty }$. Then a scaling law for steady ground effect can be deduced as
\begin{equation} C^{steady}_L = 2{\rm \pi}\beta\theta - c_1\frac{{\rm \pi}\beta^2\theta^2 D^*}{b^2+4{D^*}^2}\quad \text{where}\ c_{1} = 2.76. \end{equation}
Using the static foil numerical and experimental data from figure 8(a), the coefficient is determined to be 
$c_1 = 2.76$ by linear regression. Figure 11(a) then compares the scaling law prediction with the data, which shows that it can predict the steady ground effect data to within 
$10\,\%$ for the numerical data and 
$15\,\%$ for the experimental data.
Scaling laws for steady and unsteady ground effects show a good collapse across a wide range of 
$St$, 
$k$ and 
$D^*$. (a) The lift for a steady foil in ground effect (data from figure 8(a) with additional numerical case at 
$D^*=0.25$) is compared with the steady foil scaling law. (b) The quasi-steady lift of an unsteady pitching foil (data from figure 9 compared with its scaling relation for 
$0.25 \leq D^* \leq 1.5$, 
$0.25 \leq St \leq 0.55$ and 
$0.55 \leq k \leq 1$. The marker colour from dark blue to light yellow represents increasing 
$St$ while the marker size from small to large represents increasing 
$k$. (c) Schematic to illustrate the simple flow model for calculating the wake circulation amplification factor. (d) The wake-induced lift of an unsteady pitching foil compared with its scaling relation. Markers outlined in red represent data with 
$St \ge 0.45$. (e) Stable equilibrium altitudes compared with the scaling relation. Dashed lines present 
$20\,\%$ margins of error. These include 
$2$-D numerical data with 
$0.25 \leq St \leq 0.55$ and 
$0.55 \leq k \leq 1$, and 
$2$-D experimental data with 
$0.2 \leq St \leq 0.4$ and 
$0.6 \leq k \leq 0.9$, as well as 
$2$-D experimental data from Kurt et al. (Reference Kurt, Cochran-Carney, Zhong, Mivehchi, Quinn and Moored2019) with 
$0.2 \leq St \leq 0.6$ and 
$0.5 \leq k \leq 1$. Data points in the blue shaded area cover the Strouhal number range of biological locomotion (Taylor, Nudds & Thomas Reference Taylor, Nudds and Thomas2003; Saadat et al. Reference Saadat, Fish, Domel, Di Santo, Lauder and Haj-Hariri2017).
Returning to the time-varying formula (3.11), we can also consider an unsteady ground effect case where the foil undergoes sinusoidal pitching. If this scaling relation is then time averaged only a contribution from the second term will remain as
\begin{equation} \bar{C}_{L_{q}} =-c_{1}\frac{\beta^2\overline{{\varGamma^*}^2_{0,\infty}}D^*}{{\rm \pi}(b^2+4{D^*}^2)}. \end{equation}
Considering that 
$\theta = \theta _0 \sin (2{\rm \pi} f t)$ and 
$\theta _{0} = A^*/2=St/(2k)$ we have
$$\begin{gather} \bar{C}_{L_{q}} = c_{1}\phi_{1}, \end{gather}$$
$$\begin{gather}\phi_{1} =-\frac{{\rm \pi}^3}{8}\left(\frac{\beta^2 D^*}{b^2+4{D^*}^2}\right) \left(\frac{9}{4}+ \frac{1}{{\rm \pi}^2k^2} \right) St^2. \end{gather}$$
The coefficient 
$c_{1}=2.76$ determined from the static foil data is used for scaling the quasi-steady lift of both the steady and unsteady foil cases. Figure 11(b) presents an excellent collapse of the quasi-steady lift data over a wide range of 
$St$, 
$k$ and 
$D^*$.
3.6.2. Wake-induced lift
The wake-induced lift can be calculated as
$$\begin{gather} C_{L_{w}} \propto U^* \gamma \varGamma_{w,\infty}^*, \end{gather}$$
$$\begin{gather}U^* = 1 + U^*_{i,1} + U^*_{i,2}, \end{gather}$$
where 
$\varGamma _{w,\infty }^*$ is the normalized infinite domain wake circulation, and the effective flow velocity 
$U^*$ is composed of the velocity 
$U_{i,1}^*$, induced by an image bound vortex (representing the total bound vorticity, not the quasi-steady bound vorticity), and the velocity induced by an image wake vortex, 
$U_{i,2}^*$. To determine a scaling law for the wake-induced lift, consider a simplified flow model where a wake vortex is placed along the chord line at a distance of 
$U_\infty /(4f)$ downstream of the trailing edge (figure 11c), which gives,
\begin{align} U^*_{i,1} = \frac{-\varGamma_{\infty}^*\gamma}{\rm \pi}\left(\frac{D^*}{b^2 + 4{D^*}^2}\right)\quad \text{and}\quad U^*_{i,2} = \frac{-\varGamma_{w,\infty}^*\gamma}{\rm \pi}\left[\frac{16k^2D^*}{(1+4k)^2 + 64k^2{D^*}^2}\right], \end{align}
where 
$\varGamma ^*_{\infty }$ is the normalized infinite domain total circulation, and 
$\gamma$ is the wake circulation amplification factor. By enforcing the Kutta condition at the foil's trailing edge, the wake circulation amplification factor is calculated as
\begin{equation} \gamma = \frac{1+4k(1-b)}{1+4k(1-b) - \dfrac{(1-b)^2}{(1-b)^2+4{D^*}^2} - \dfrac{4k(1-b)}{1+64k^2{D^*}^2}}. \end{equation}Following Moored & Quinn (Reference Moored and Quinn2018), we have
$$\begin{gather} \varGamma_{w,\infty}^* \propto \varGamma_{0,\infty}^* \frac{k^*}{k^*+1} \quad \text{and}\quad k^* = \frac{k}{1+ 4St^2}, \end{gather}$$
$$\begin{gather}\varGamma_{\infty}^* \propto \varGamma_{0,\infty}^*. \end{gather}$$
Similar to the quasi-steady lift, after time averaging the wake-induced lift, only terms associated with 
${\varGamma ^*}^2_{0,\infty }$ remain, which gives
$$\begin{gather} \bar{C}_{L_{w}} = c_{2}\phi_{2} + c_{3}\phi_{3}, \end{gather}$$
$$\begin{gather}\phi_{2} =-\gamma^2St^2 \left(\frac{9}{4}+ \frac{1}{{\rm \pi}^2k^2} \right)\left(\frac{D^*}{b^2+ 4{D^*}^2}\right)\left(\frac{k^*}{k^*+1}\right), \end{gather}$$
$$\begin{gather}\phi_3 =-\gamma^2St^2 \left(\frac{9}{4}+ \frac{1}{{\rm \pi}^2k^2} \right)\left[\frac{16k^2D^*}{(1+4k)^2 + 64k^2{D^*}^2}\right]\left(\frac{k^*}{k^*+1}\right)^2. \end{gather}$$
By using the wake-induced lift data from figure 4, the coefficients were determined by linear regression to be 
$c_{2} = -33.00$ and 
$c_{3} = 83.35$. Figure 11(d) shows that a good collapse of the data can be achieved with the exception of data with 
$St \ge 0.45$ (outlined in red). As shown in figure 7(a), the wake-induced lift is significantly altered by a wake instability effect that occurs when 
$St \ge 0.45$ even for out-of-ground effect foils. Since this wake instability effect is not accounted for in the scaling laws, the data points with 
$St \ge 0.45$ (figure 11d) begin to show deviation from the scaling law prediction. Details of the scaling laws derivation, such as calculation of circulation amplification factors, velocities induced by the image system and time-averaged circulation square, can be found in Appendix A.
3.6.3. Total lift
Now, the quasi-steady lift and wake-induced lift can be summed to acquire the total lift, and the total lift becomes a function of 
$St$, 
$k$ and 
$D^*$
\begin{equation} \bar{C}_{L} = c_{1}\phi_{1} + c_{2}\phi_{2} + c_{3}\phi_{3}= f \left(St, k, D^*\right). \end{equation}
Then, for a given pair of 
$St$ and 
$k$, the equilibrium altitude, 
$D^*_{eq}$, can be determined by setting (3.25) to zero, which gives
\begin{align} & c_{1}\frac{{\rm \pi}^3}{8}\left(\frac{\beta^2 }{b^2+4{D^*_{eq}}^2}\right) + c_{2}\gamma^2 \left(\frac{1}{b^2+ 4{D^*_{eq}}^2}\right)\left(\frac{k^*}{k^*+1}\right) \nonumber\\ &\quad +c_{3} \gamma^2 \left[\frac{16k^2}{(1+4k)^2 + 64k^2{D^*_{eq}}^2}\right]\left(\frac{k^*}{k^*+1}\right)^2=0. \end{align}
The equilibrium altitude cannot be written in a closed-form explicit solution, but rather (3.26) must be solved implicitly. In this way, (3.26) represents an implicit scaling law for the equilibrium altitude of a pitching foil in ground effect. Figure 11(e) presents the scaling law prediction for the equilibrium altitude. The scaling law can predict the equilibrium altitude to within 
$20\,\%$ error except for the data with 
$St \ge 0.45$. Moreover, the predicted equilibrium altitude for 
$St\ge 0.45$ is shown to decrease with increasing 
$St$. The deviation in the prediction of equilibrium altitudes at 
$St \ge 0.45$ is due to the deviation in the wake-induced lift prediction, which is caused by neglecting the wake instability effect in the wake-induced lift scaling law. At 
$St \ge 0.45$ the foil has a non-zero time-averaged wake-induced lift as 
$D^*$ goes to infinity. However, the time-averaged wake-induced lift always goes to zero in the scaling law as 
$D^*$ goes to infinity. Data points in the blue shaded area in figure 11(e) cover the Strouhal number range of biological locomotion that is 
$0.2\leq St \leq 0.4$ (Taylor et al. Reference Taylor, Nudds and Thomas2003; Saadat et al. Reference Saadat, Fish, Domel, Di Santo, Lauder and Haj-Hariri2017), so it can be observed that the scaling law captures the trend of the biologically relevant data. For the numerical data with 
$St<0.45$ a 
$5\,\%$ and 
$10\,\%$ error in predicting the quasi-steady and wake-induced lift, respectively, leads to a 
$20\,\%$ error in predicting the equilibrium altitude. This indicates that the equilibrium altitude relies on a delicate balance and is sensitive to small changes in the lift components. This is also observed by Liu et al. (Reference Liu, Zhong, Han, Moored and Quinn2023), who found that a pitch bias angle equal to 
$7\,\%$ of the amplitude can shift the equilibrium altitude of a pitching foil by 
$28\,\%$.
4. Conclusions
In this article we decomposed the lift of a near-ground pitching hydrofoil into its added mass, wake-induced and quasi-steady components, and determined that their time-averaged values are always zero, positive and negative, respectively, across all ground proximities. This shows that both stable and unstable equilibrium altitudes are generated by a balance between positive wake-induced lift and negative quasi-steady lift while the added mass lift does not play a role.
Using both simulations and experiments, detailed analyses are provided to illustrate the three lift components’ near-ground behaviour. In ground effect the negative peak of the added mass lift is amplified more than the positive peak, however, the negative peak also becomes sharper while the positive peak becomes broader compared with a foil out of ground effect. This leads to a scallop theorem result, where the added mass lift precisely integrates to zero over one oscillation cycle. Through a control volume analysis, the wake-induced lift is discovered to be positive due to the positive lift from the pressure difference on the control volume outweighing the negative lift from the net momentum flux. This leads to the upwards deflected wake behind a near-ground swimmer being responsible for a positive wake-induced lift. The quasi-steady lift is further analysed through the lens of the Kutta–Jowkouski theorem. It is observed that the velocity induced by the positive (negative) quasi-steady circulation bounded to the ground image foil enhances (reduces) the effective flow velocity at the foil's leading edge. For a pitching foil, this amplifies the negative quasi-steady lift during the upstroke and reduces the positive quasi-steady lift during the downstroke leading to a net negative quasi-steady lift. Additionally, the lift decomposition was applied to 3-D pitching foil simulations of varying aspect ratio. It was demonstrated that equilibrium altitudes of 3-D hydrofoils are also generated by a balance between positive wake-induced lift and negative quasi-steady lift. Moreover, as the aspect ratio decreases, the magnitude of the wake-induced lift drops off faster than the magnitude of the quasi-steady lift, leading to the disappearance of an equilibrium altitude when 
${A{\kern-4pt}R} < 1.5$.
Using these insights, scaling laws of the quasi-steady lift, wake-induced lift and equilibrium altitudes are developed. Additionally, a simple scaling law for the lift of a steady foil in ground effect is discovered. The scaling laws show good agreement to both numerical and experimental data by predicting the equilibrium altitudes to within 
$20\,\%$ of their values as long as 
$St<0.45$. For 
$St \ge 0.45$, significant effects from wake instability, not accounted for in the scaling laws, arise that alter the wake-induced lift. These results not only provide key physical insights and scaling laws for unsteady ground effect, but also for two schooling hydrofoils in a side-by-side formation with an out-of-phase synchronization.
Funding
This work was supported by the National Science Foundation under Program Director Dr R. Joslin in Fluid Dynamics, award number 1921809 as well as by the Office of Naval Research under Program Director Dr R. Brizzolara on MURI grant number N00014-22-1-2616.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of scaling laws
A.1. Quasi-steady circulation amplification factor
When the foil is near the ground (figure 12a), in order for the Kutta condition to be preserved at the foil's trailing edge, an additional bound circulation of the real foil must arise to cancel the transverse velocity induced by the image foil, which gives
\begin{equation} (\beta-1)V_{1} + V_{2}=0, \end{equation}with
\begin{equation} V_{1}=\frac{\varGamma_{0,\infty}}{2{\rm \pi}(1-b)c}\quad \text{and}\quad V_{2} =-\frac{\varGamma_{0,\infty}\beta}{2{\rm \pi}}\left[\frac{(1-b)c}{(1-b)^2c^2+4D^2}\right], \end{equation}
where 
$V_{1}$ is the transverse velocity induced by the bound vortex of the real foil out of ground effect, and 
$V_{2}$ is the transverse velocity induced by the bound vortex of the image foil.
Schematics for calculating (a) the quasi-steady circulation amplification factor, 
$\beta$, and (b) the wake circulation amplification factor, 
$\gamma$.
Then, 
${\varGamma _{0,\infty }}/{2{\rm \pi} c}$ can be divided from both sides, which gives
\begin{equation} \frac{\beta - 1}{1-b} - \beta\left[\frac{(1-b)}{(1-b)^2+4{D^*}^2}\right] = 0. \end{equation}After organization we have
\begin{equation} \beta = 1 + \frac{(1-b)^2}{4{D^*}^2}. \end{equation}A.2. Wake circulation amplification factor
When the foil is near the ground (figure 12b), in order for the Kutta condition to be preserved at the foil's trailing edge, additional bound and wake circulations of the real foil must arise to cancel transverse velocities induced by the image foil, which gives
\begin{equation} (\gamma-1)V_{3} + (\gamma-1)V_{4} + V_{5} + V_{6} = 0, \end{equation}with
$$\begin{gather} {V_3} = \frac{\varGamma_{w,\infty}}{2{\rm \pi}(1-b)c},\quad V_{4} = \frac{2\varGamma_{w,\infty}\,f}{{\rm \pi} U_{\infty}}, \end{gather}$$
$$\begin{gather}V_{5} =-\frac{\varGamma_{w,\infty}\gamma(1-b)c}{2{\rm \pi}[(1-b)^2c^2+4D^2]},\quad \text{and}\quad V_{6} =-\frac{2\varGamma_{w,\infty}\gamma f U_{\infty}}{{\rm \pi}(U^2_{\infty}+64D^2f^2)}, \end{gather}$$
where 
$V_{3}$ and 
$V_{4}$ are transverse velocities induced by the bound and wake vortex of the real foil out of ground effect, and 
$V_{5}$ and 
$V_{6}$ are transverse velocities induced by the bound and wake vortex of the image foil.
Then 
${\varGamma _{w,\infty }}/{2{\rm \pi} }$ can be divided from both sides, which gives
\begin{equation} \frac{\gamma- 1}{(1-b)c} + \frac{4f(\gamma-1)}{ U_{\infty}} - \frac{\gamma(1-b)c}{(1-b)^2c^2+4D^2} - \frac{4\gamma f U_{\infty}}{U^2_{\infty}+64D^2f^2}=0. \end{equation}After organization we have
\begin{equation} \gamma\left[ \frac{1}{1-b} + 4k - \frac{1-b}{(1-b)^2+ 4 {D^*}^2}- \frac{4k}{1+64{D^*}^2k^2} \right] = \frac{1}{1-b} + 4k. \end{equation}
Then 
$\gamma$ can be calculated as
\begin{equation} \gamma = \frac{1+4k(1-b)}{1+4k(1-b) - \dfrac{(1-b)^2}{(1-b)^2+4{D^*}^2} - \dfrac{4k(1-b)}{1+64k^2{D^*}^2}}. \end{equation}A.3. Calculation of time-averaged circulation square
The normalized quasi-steady circulation is
\begin{equation} \varGamma_{0,\infty}^* = {\rm \pi}\theta +\frac{3{\rm \pi} c \dot{\theta}}{4U_{\infty}} = {\rm \pi}\theta_{0}\sin(2{\rm \pi} f t) + \frac{3{\rm \pi}^2 c f \theta_{0 }\cos(2{\rm \pi} f t)}{2U_{\infty}}, \end{equation}
which, considering 
$\theta _{0} = A^*/2=St/(2k)$, becomes
\begin{equation} \varGamma_{0,\infty}^* = \frac{{\rm \pi} St}{2k}\sin(2{\rm \pi} f t) + \frac{3{\rm \pi}^2St}{4}\cos(2{\rm \pi} f t). \end{equation}
The time-averaged square of the normalized quasi-steady circulation, 
$\overline {{\varGamma ^*}^2_{0,\infty }}$, is calculated as
\begin{equation} \overline{{\varGamma^*}^2_{0,\infty}} = f\int_{0}^{1/f} \left[\frac{{\rm \pi} St}{2k}\sin(2{\rm \pi} f t) + \frac{3{\rm \pi}^2St}{4}\cos(2{\rm \pi} f t)\right]^2 \text{d}t. \end{equation}After expansion we have
\begin{align} \overline{{\varGamma^*}^2_{0,\infty}} = f\int_{0}^{1/f} \frac{{\rm \pi}^2 St^2}{4k^2}\sin^2(2{\rm \pi} f t) + \frac{3{\rm \pi}^3St^2}{4k}\sin(2{\rm \pi} f t)\cos(2{\rm \pi} f t) + \frac{9{\rm \pi}^4St^2}{16}\cos^2(2{\rm \pi} f t) \,\text{d}t. \end{align}The second term in the above equation integrates to zero. Then we have
\begin{equation} \overline{{\varGamma^*}^2_{0,\infty}} = f\int_{0}^{1/f} \frac{{\rm \pi}^2 St^2}{4k^2}\sin^2(2{\rm \pi} f t) + \frac{9{\rm \pi}^4St^2}{16}\cos^2(2{\rm \pi} f t)\,\text{d}t = \frac{{\rm \pi}^4St^2}{8}\left(\frac{9}{4} + \frac{1}{{\rm \pi}^2k^2}\right). \end{equation}A.4. Calculation of the induced flow velocity
The horizontal velocity induced by an image bound circulation, 
$\varGamma$, at the foil's leading edge is calculated as
\begin{equation} U_{i,1} = \frac{\varGamma}{2{\rm \pi} \sqrt{b^2c^2 + 4D^2}} \left(\frac{2D}{\sqrt{b^2c^2+4D^2}}\right) = \frac{\varGamma D}{{\rm \pi}(b^2c^2+4D^2)}. \end{equation}After normalization we have
\begin{equation} U^*_{i,1} = \frac{\varGamma^*}{\rm \pi} \left(\frac{D^*}{b^2+4{D^*}^2}\right). \end{equation}
The horizontal velocity induced by a image wake circulation, 
$\varGamma$, at the foil's leading edge is calculated as
\begin{align} U_{i,2} &= \frac{\varGamma}{2{\rm \pi} \sqrt{4D^2+ \left(c + \dfrac{U_{\infty}}{4f}\right)^2}} \left(\frac{2D}{\sqrt{4D^2+ \left(c + \dfrac{U_{\infty}}{4f}\right)^2}}\right) \nonumber\\ &=\frac{\varGamma D}{{\rm \pi}\left[4D^2+ \left(c + \dfrac{U_{\infty}}{4f}\right)^2\right]}. \end{align}After normalization we have
\begin{equation} U^*_{i,2} = \frac{\varGamma^*}{\rm \pi}\left[\frac{16k^2D^*}{(1+4k)^2 + 64k^2{D^*}^2}\right]. \end{equation}
Note that 
$\varGamma$ is a general representation for vortices at locations of the image bound or wake vortex, and it is replaced by circulations associated with 
$\varGamma _{0,\infty }$ or 
$\varGamma _{w,\infty }$ in the article accordingly.
