The Sensitivity of the Fitch Wind Farm Parameterization to a Three-Dimensional Planetary Boundary Layer Scheme

. Wind plant wake impacts can be estimated with a number of simulation methodologies, each with its own ﬁdelity and sensitivity to model inputs. In turbine-free mesoscale simulations, hub-height wind speeds can signiﬁcantly vary with the choice of a planetary boundary layer (PBL) scheme. However, the sensitivity of wind plant wakes to a PBL scheme has not been explored because, until now, wake parameterizations were only compatible with one PBL scheme. We couple the Fitch wind farm parameterization with the new NCAR 3DPBL scheme and compare the resulting wakes to those simulated with a 5 widely used PBL scheme. First, we simulate a wind plant in a pseudo-steady state under idealized stable, neutral, and unstable conditions with two PBL schemes: MYNN and the NCAR 3DPBL. For these idealized scenarios, MYNN consistently predicts internal wakes that are 0.25–1.5 m s − 1 stronger than internal 3DPBL wakes. However, because MYNN predicts stronger inﬂow winds than 3DPBL, MYNN predicts average capacity factors that are as large as 13 percentage points higher than with the 3DPBL, depending on the stability. To extend this sensitivity study, we conduct a month-long case study with both PBL 10 schemes centered on the Vineyard Wind 1 lease area in the mid-Atlantic United States. Under stable and unstable conditions averaged across the month, MYNN again predicts stronger internal waking—by about 0.25 m s − 1 . However, again due to stronger plant inﬂow wind speeds in MYNN, the 3DPBL generates 4.7%–7.8% on the turbine build-out scenario. Differences between PBL schemes can be even larger for individual instances in time. These simulations suggest that PBL schemes represent a meaningful source of modeled wind resource uncertainty; therefore, we recommend incorporating PBL variability into future wind plant planning sensitivity studies as well as wind forecasting studies. time-averaged visualizations. MYNN shows a greater prevalence of winds at rated power (Region 510 III) under stable and unstable conditions. In neutral conditions, both MYNN and 3DPBL winds are predominantly in Region II; therefore, we expect less wind resource (larger wake effects on power production) in this stability condition, at least for the month of

turbines, finding that the bias shifted between 4 MW and 10 MW when vertical grid resolution was varied and 1 MW and 7 MW when the forcing product was varied. 90 Mangara et al. (2019) simulated three weeks of onshore winds in flat terrain in Changyi, China. They examined the impact of horizontal and vertical resolution on modeled wakes. They calculated relative wind speed deficit in the most-downwind grid cells that contained turbines, and they found that 19-day averaged hub-height relative wind speed deficits ranged between approximately 4% and 10% when horizontal resolution was varied and approximately 5% and 9% when vertical resolution was varied. They also found that added averaged TKE values could significantly vary with grid resolution, finding that added 95 hub-height TKE in the last row of the plant varied between approximately 0.15-0.50 m 2 s −2 , depending on the grid resolution. Shepherd et al. (2020) conducted a yearlong simulation over Iowa with the Fitch WFP and the EWP. They calculated vertical profiles of median velocity deficits inside the Pomeroy wind plant cluster, finding that Fitch could predict winds that were 1 m s −1 stronger than in the EWP. The two schemes also predicted different heights for maximum wind speed deficit, where the EWP produced the strongest deficits within the rotor disk, whereas, at times, Fitch produced wake maxima that lofted 100 above the rotor disk. External wakes (here defined as the distance at which wind speed deficits at all heights were smaller than 0.2 m s −1 ) have a similar length for both PBL schemes, although wind speed deficits inside the external wake were stronger with Fitch because of its stronger internal waking. The wind speed deficits correspondingly produce different power losses, where the EWP predicts 5% more total annual power production than Fitch. Finally, due to the fundamentally different manner in which TKE is treated within the WFPs, median values of TKE can be 1 m 2 s −2 stronger with Fitch than with the EWP inside 105 the plant. Pryor et al. (2020) simulated 9 months of winds over Iowa with a number of model configurations. When they increased the number of vertical levels from 41 to 57, monthly mean hub-height wind speeds increased between 0.2 m s −1 and 0.6 m s −1 in turbine-free simulations. However, when horizontal grid resolution was doubled from 4 km to 2 km in simulations with 41 vertical levels, the turbine-free hub-height wind speed distributions were statistically similar. The differences in maximum 110 hub-height TKE enhancement from the Fitch WFP in each month of the study varied between 0.2 m 2 s −2 and 0.3 m 2 s −2 , depending on the grid resolution. Median gross capacity factors simulated with the Fitch WFP were 38%-41%, depending on While the values of empirical constants are different, MYNN and the quasi-3DPBL use the same formulation to parameterize turbulent momentum, heat, and moisture fluxes. For example, they parameterize the vertical flux of the u-component of wind speed as where L is the master length scale, q is √ 2 TKE, S m is a stability function, and U is zonal velocity (Mellor and Yamada,160 1982).
The two PBL schemes also use identical formulations for components of the TKE budget. We briefly recount their formulations here, because we discuss TKE budgets in Sec 3.6. Precise details of each formulation, as well as definitions of each variable, can be found in Nakanishi and Niino (2009), Olson (2019), and ?. Both PBL schemes calculate shear-generated TKE as P s = − uw ∂U ∂z + vw ∂V ∂z , buoyancy-generated TKE as P b = −βg i wθ v , the diffusion of TKE as ∂ ∂x k LqS q ∂q 2 ∂x k , and 165 TKE dissipation as ε = q 3 /B 1 L. While the formulations are identical, the budget components take on different values in each PBL scheme due to the value of different empirical constants and the numerical implementation of each calculation.

Integration of the Fitch WFP with the 3DPBL
To simulate wakes with the 3DPBL, we first integrated the Fitch WFP with the 3DPBL inside the WRF code. The Fitch WFP modifies flow in two key manners (Fitch et al., 2012;Archer et al., 2020): by slowing winds and by adding TKE In the above equations, k is the vertical level that intersects the rotor, A k is the area of the rotor on this vertical level, C T is 175 the thrust coefficient, U k is the wind speed, u k is the zonal wind, v k is the meridional wind, and z k is the height. The turbulence coefficient C T KE is calculated as the difference between the thrust coefficient C T and the power coefficient C P . The thrust and power coefficients are functions of wind speed that are unique to a particular wind turbine, and their values are specified in the input file wind-turbine.tbl. The coefficient c a was introduced by Archer et al. (2020) to empirically modify the amount of explicit TKE addition and, in this study, we either set it to 0 or 1. The major challenge in integrating the Fitch WFP and the 3DPBL is that the 3DPBL code is housed in the dynamics (dyn_em/ ) part of the code, as opposed to the physics (phys/ ) part of the code where most other PBL schemes reside. As such, the codebase was modified to account for the user-selected PBL scheme. A call to the Fitch WFP's dragforce subroutine was added to the end of dyn_em/module_first_rk_step_part2.F. When called for the 3DPBL, the velocity tendencies and TKE tendencies are additionally scaled by the column-mass in order to match the identical scaling that happens to the phys/ -185 calculated tendencies earlier within dyn_em/module_first_rk_step_part2.F. Additionally, whereas the Fitch WFP code modifies the MYNN TKE field directly (including a timestep factor of ∂t), the new code modifies the 3DPBL TKE tendency field (omitting a factor of ∂t and letting the rest of the code carry out the time integration).

Configuration of Idealized Simulations
First, we carry out a series of idealized simulations to study the effect of the PBL scheme on simulated wake dynamics in 190 a simple flat terrain environment. All simulation inputs can be found on Zenodo (https://doi.org/10.5281/zenodo.5565399).
We use the neutral idealized simulations of Fitch et al. (2012) as inspiration for our simulations, but we make a number of modifications. All simulations use two domains, each 202-by-202 grid points in the horizontal. MYNN is always used in the outer domain, whereas the inner domain is either MYNN or the 3DPBL. The outer domain uses a horizontal grid spacing of 3 km and a timestep of 9 seconds, whereas the inner domain uses a horizontal grid spacing of 1 km and a timestep of 3 seconds.

195
The vertical grid uses 81 cells, up to a height of 20 km. Vertical grid stretching is employed to provide finer resolution near the surface, thereby allowing 28 vertical levels below a height of 300 m, following the recommendation of Tomaszewski and Lundquist (2020) (Rosencrans et al., in prep), and are smaller than typical values over land. Simulations for each stability case are initialized with a neutral temperature profile of 285 K within the boundary layer up to 500 m. The boundary layer is capped with a two-layer inversion: a strong inversion (5 K warming between 500 m and 600 m) and a weaker inversion (3 K km −1 lapse rate above 600 m). One turbine-free simulation is 205 spun up for three days for each stability and each PBL scheme, after which, pseduo-steady conditions are achieved. This long spinup is done to dampen inertial oscillations, as in Fitch et al. (2012). After spin up, the boundary-layer height is approximately 600 m in the stable simulations, 750 m in the neutral simulations, and 1,150 m in the unstable simulations.
After spinning up turbine-free simulations, we run three cases of simulations for each of the stabilities and each of the PBL schemes (Table 1). The first case is simply a continuation of the turbine-free simulations and is referred to as the no-wind-farm 210 (NWF) case. The second case (100TKE) includes a 10-by-10 grid of 12-MW International Energy Agency (IEA, Beiter et al., 2020) reference offshore wind turbines placed in the center of the inner domain. The turbine hub-height is 138 m, and the rotor diameter is 215 m. Cut-in speed is 3 m s −1 , rated speed is 10.9 m s −1 , and cut-out speed is 30 m s −1 . Turbines are placed 2 km apart, which is close to the 1 nautical mile spacing used in the real simulations. In this case, 100% of explicit TKE is generated by the Fitch WFP (c a = 1). In the third case (0TKE), we explore the sensitivity to explicitly added TKE by duplicating the 215 setup of the second case, but turning off explicit TKE generation (c a = 0).  Next, to ground this analysis in a real-world application, we study wake sensitivity centered on planned offshore wind plants off the U.S. east coast (Fig. 1) both of the PBL schemes ( Table 2). The first case simulates a turbine-free atmosphere and is referred to as the NWF case. We focus on Vineyard Wind because, at the time the simulations for this study were initiated, it was likely to be first 100+ MW project in US offshore water.

Configuration of Real
We simulate winds for the month of August 2020. We chose this period because of its high electricity demand (Livingston 235 and Lundquist, 2020). We start the simulations on July 30, 2020, to allow for 48 hours of spin-up, and we omit this data from analysis. The domain used in this study is identical to the one used by the 20-year NREL mid-Atlantic analysis available at NREL (2020) as well as Rosencrans et al. (in prep). The outer domain at 6-km grid spacing has 196 grid points in the west-east direction and 122 grid points in the south-north direction and uses a timestep of 18 seconds. The inner domain of 2-km grid spacing has 466 grid points in the west-east direction and 259 in the south-north direction and uses a timestep of 6 seconds. 240 We save data from the inner domain every 5 minutes. As with the idealized simulations, the outer domain for all simulations uses the MYNN PBL scheme, whereas the inner domain either uses MYNN or the 3DPBL. TKE advection is turned on for all MYNN domains. All domains use Thompson microphysics (Thompson et al., 2008), RRTMG for radiation (Iacono et al., 2008), Jiménez-modified Monin-Obukhov for the surface layer scheme (Jiménez et al., 2012), and the unified Noah landsurface model (Tewari et al., 2004). Horizontal turbulent mixing is carried out with a Smagorinsky-diffusion style approach in 245 MYNN, whereas it is calculated from the horizontal turbulent flux divergence in the 3DPBL. ERA5 (Hersbach et al., 2020) provides atmospheric forcing, and OSTIA at 1/10 • resolution provides sea surface temperatures (Donlon et al., 2012).
We select turbines and turbine spacing that are consistent with current expected standards of offshore U.S. wind plants. As was done in the idealized simulations, we use the 12-MW IEA reference turbine in the real simulations. This is similar to the 62 13-MW turbines that are slated for operation in Vineyard Wind 1 (Vineyard Wind, 2021). All modeled lease areas are 250 fully built out in order to study the maximum possible power production and wake strength. The turbine spacing is identical to that used by Rosencrans et al. (in prep). They are spaced 1 nautical mile apart, which is the same spacing that will be used in Vineyard Wind 1. In total, 177 turbines are modeled in the VW-ONLY simulations, and 1,418 turbines are modeled in the LEASE simulations. We note that the official layout of the Vineyard Wind 1 site was announced after our simulations were completed. The official layout will only include turbines in the northern half of Vineyard Wind 1. As such, our study will 255 overestimate the magnitude of internal waking at Vineyard Wind 1 in the years immediately after it is built.

Turbine-Free Conditions
Before comparing the impacts of the turbines on the wind flow, it is necessary to characterize the baseline wind speeds and TKE in the turbine-free simulations. All simulations are forced with 10 m s −1 geostrophic winds, but wind speed profiles are Here, MYNN hub-height wind speeds are 0.22 m s −1 (3%) stronger than in the 3DPBL.
In these neutral and stable conditions, the two PBL schemes can produce different amounts of TKE in NWF simulations ( Fig. 2d,  The differences and similarities in the wind profiles of the NWF simulations will dictate the comparison of the wake effects 285 between the PBL schemes in turbine-including simulations. Differences in wake magnitude will largely be governed by differences in the wind speeds entering the plant as differences in wake recovery processes, which are linked to parameterizations of turbulent fluxes (Gupta and Baidya Roy, 2021). The turbine thrust coefficient is nearly identical for all simulations based on hub-height wind speeds, 0.803-0.806, simplifying the comparison. Additionally, the two PBL schemes produce similar NWF wind speed profiles under neutral and stable conditions. Thus, in these two conditions, we expect that wake differences 290 will primarily arise from differences in the way that turbulent fluxes are modeled within each PBL scheme. Under unstable conditions, however, the differences in NWF wind speeds will additionally play a major role in wake dynamics. Finally, we note that NWF hub-height wind speeds increase in strength when moving from stable to neutral to unstable conditions. This wind speed variability will impact power production and losses when comparing between stability conditions.

295
Internal wakes are sensitive to the choice of PBL scheme and stability (Fig. 3). We quantify internal wake strength by finding the average hub-height wind speeds within the plant in the turbine-including simulations ("WFP," which is a generic stand-in for "100TKE" or "0TKE") relative to hub-height winds in the turbine-free simulations ("NWF"). We also calculate the percentage of wind speed loss with reference to the NWF winds inside the plant, and we calculate percentage point [pp] differences as differences between percentage values. The idealized MYNN simulations consistently show internal wakes that are stronger 300 than their 3DPBL counterparts at hub-height. These differences primarily arise from two sources-differences in inflow wind speed and differences in turbulent mixing formulations. The neutral and stable simulations have similar inflow wind speeds.
As such, for neutral and stable cases, the differences in internal wakes are more driven by differing wake recovery rates than by differing inflow wind speed. Recovery rates differ between PBL schemes because different empirical constants are used to calculate turbulent fluxes and because the PBL schemes have different ambient vertical profiles of TKE (particularly in stable conditions. Additionally, we note that in these idealized simulations, internal wakes are typically more sensitive to PBL scheme choice than the presence or absence of explicitly added TKE. While external wakes are more nebulous to characterize, they also show sensitivity to atmospheric stability. There is no 320 singular standard approach that is used to characterize external wakes (Fischereit et al., 2021), so we adopt two approaches: by (an effect which is observed in these simulations and will be discussed momentarily). Thus, the larger TKE under neutral MYNN also helps erode the neutral MYNN wakes in these simulations.
We additionally characterize external wakes in the idealized simulations by using the e-folding distance. Internal wakes in the 3DPBL have values closer to 0.5 m s −1 ; thus, the 0.5 m s −1 deficit contour is less useful to characterize external wakes. Instead, we calculate the e-folding contour as 1/e times the average internal wake strength, or about 36% (Fitch et al.,335 2012). This metric characterizes external wakes relative to internal wake magnitude. In the 100TKE 3DPBL simulations, the stable and neutral e-folding distance extends beyond 80 km, whereas it is 25 km in the unstable simulations. Thus, the unstable 100TKE 3DPBL simulations show quicker wake recovery, as expected. By contrast, the e-folding distance in the 0TKE 3DPBL simulations show little sensitivity to stability, as all the wake lengths are 12 km-17 km, perhaps due to similar values of ambient turbulence in the 3DPBL simulations ( Fig. 2d-f).

340
Depending on the presence of TKE addition and stability, MYNN and the 3DPBL can produce similar external wakes.
The addition of explicit TKE inside the plant consistently either has negligible impact on external wakes or it extends the external wakes further. This is somewhat counterintuitive on first glance, as the addition of explicit TKE weakens internal wakes. However, the external wake length growth may be associated with a slight decrease in turbulent mixing due to a small decrease in TKE downwind of the plant (as observed in Fig. 6g and Fitch et al. (2012) We briefly digress from the discussion on wakes to note that upwind blockage (Schneemann et al., 2021

Vertical Structure of Wind Speed Deficits
While hub-height winds are particularly important to quantify, it is also helpful to characterize wakes over the vertical extent of the rotor disk. We calculate the wind speed deficit averaged along the y-extent (predominantly crosswind) of the plant for 370 each simulation (Fig. 4).
Just as MYNN consistently predicted stronger hub-height wind speed deficits, it also tends to predict stronger internal wakes over the extent of the rotor disk. MYNN tends to predict maximum wind speed deficits inside the plant that are 0.25-0.50 m s −1 stronger than those of the 3DPBL. The one exception to this is the stable 100TKE simulations (Fig. 4a,

Difference in Momentum Tendencies
In large part, the two PBL schemes produce different wind speed deficits in their wakes because the schemes parameterize turbulent fluxes differently, as we illustrate here. The u and v components of wind speed are modified by mechanisms such as advection of the mean wind, the Coriolis force, and the divergence of the turbulent momentum fluxes (Stull, 1988, Eqn. 5.1a  We expect all these terms, aside from the divergence of turbulent momentum fluxes, to be similar for both MYNN and the 3DPBL, as the NWF wind speeds are similar but two PBL schemes parameterize momentum fluxes uniquely. We calculate the u-tendency due to the turbulent flux divergence as the vertical derivative of u w , neglecting the horizontal components of flux divergence because they are significantly smaller in the 3DPBL than u w , and they are not computed in the MYNN parameterization. We also omit visualizations of v-tendency because they are substantially smaller than the u-tendency in these 395 idealized simulations forced with a u geostrophic wind. We investigate the relationship of wind speeds and turbulent fluxes between the two PBL schemes in the neutral 100TKE simulations by comparing two fields in the wakes-the wind speed deficits and the turbulent flux divergence u-tendency "deficits" (Fig. 5). The u-tendency deficits are defined as tendencies in the turbine-free simulations subtracted from tendencies in the turbine-including simulations.
The differences in tendency deficits between the two PBL schemes drive the differences in the wind speed wakes. As winds 400 advect primarily along the x-direction, wind speed magnitudes are modified by the tendency. For example, the u-tendency is more negative for MYNN in the first few columns of the plant. Correspondingly, the MYNN wind speed deficits are more negative. They become more negative as the winds pass over regions where MYNN continues to have a more negative tendency, and they become more positive as they pass over regions where MYNN has a more positive tendency. Thus, the modeled wind speed deficits in the wake of a plant depend on how the PBL scheme parameterizes turbulent momentum fluxes.

Total TKE
While wind speed deficits are often significantly sensitive to PBL scheme and stability, shifts in TKE have a more variable response to these factors (Fig. 6). All 0TKE simulations tend to show minor TKE generation (0.1-0.3 m 2 s −2 ) at the top of the rotor disk where shear is the 420 strongest, and they also show small TKE deficits (0.1-0.2 m 2 s −2 ) below hub height relative to the NWF simulations. The one exception is the unstable 0TKE MYNN simulation, which shows negligible TKE changes when turbines are included. Thus, while TKE magnitudes in these mesoscale simulations are sensitive to the amount of TKE explicitly added by the Fitch wind farm parameterization, they are relatively insensitive to shifts in atmospheric stability and PBL scheme. Additionally, we note that TKE fields appear "smoother" in the 3DPBL simulations than they do in the MYNN simulations. This smoothing likely 425 arises from differences in the way that TKE advection, empirical constants, and horizontal mixing are handled within each PBL scheme.

TKE Budget
The changes in total TKE fields fundamentally arise from differences in how different components of the TKE budget (Stull, 1988, Eqn. 5.1a) are parameterized. We calculate profiles of four budget components over the horizontal extent of the plant 430 (Sec. 2.1): shear-generated TKE, buoyancy-generated TKE, vertical turbulent transport of TKE, and TKE dissipation. We visualize the value of the of the turbine-free profiles subtracted from the turbine-including profiles (Fig. 7).    column 1 is shear-generated TKE, column 2 is buoyancy generation or destruction of TKE, column 3 is vertical turbulent transport and column 4 is TKE dissipation. Row 1 shows these components in stable conditions, row 2 shows neutral conditions, and row 3 shows unstable conditions. Each "100TKE" and "0TKE" profile is calculated with respect to the NWF profiles as WFP-NWF. Larger differences between MYNN and the 3DPBL emerge in the 100TKE simulations, although these differences can still be quite small, again explaining similar total TKE fields in the vicinity of the plants. While shear-generated TKE is near-identical for both PBL schemes in the neutral and unstable simulations, the 3DPBL produces substantially more shear-440 generated TKE above the rotor disk than MYNN under stable conditions (Fig. 7a). This greater shear-generated TKE translates to overall larger values of TKE in this region for the 3DPBL (Fig. 6). While profiles of the vertical turbulent transport of TKE are noisy, both PBL schemes tend to have negative TKE turbulent transport in the upper half of the rotor disk and positive TKE turbulent transport in the lower half of the disk. The magnitude tends to be stronger in MYNN than in the 3DPBL, regardless of stability. The 100TKE simulations can also show larger differences in TKE dissipation, particularly in unstable conditions 445 (Fig. 7j). TKE dissipation is parameterized in both PBL schemes with a dependency on q 3 , and as such, small differences in total TKE can lead to large discrepancies between the two PBL schemes.

Power
Power production and losses due to internal waking change with PBL scheme (Fig. 8). We calculate the capacity factor for each turbine, the average capacity factor of the plant, and the average power deficit due to internal wakes with reference to the 450 NWF hub-height wind speed. Capacity factor is defined as the ratio of actual power output relative to the maximum possible power output. Average power production throughout the plant is a function of NWF hub-height wind speed and wake recovery rate (which is governed by turbulent fluxes in each PBL scheme). Under neutral and stable conditions, the 3DPBL consistently predicts larger capacity factors than MYNN by 1 pp-2 pp. This behavior occurs because the two PBL schemes have similar inflow speeds, but MYNN has stronger hub-height internal waking (Sec. 3.2). The similar NWF winds in these stabilities 455 (differing by ∼0.25 m s −1 ) allow us to isolate the effect of differing wake recovering rates between PBL schemes. Capacity factor losses due to internal waking in these stabilities are 3 pp-5 pp weaker in the 3DPBL than in MYNN. The 3DPBL power losses range from 6.6 pp-10.8 pp, whereas the MYNN power losses range from 9.9 pp-15.7 pp. Thus, in neutral and stable conditions, the 3DPBL shows a quicker wake recovery rate inside the plant.
By contrast, unstable NWF MYNN hub-height wind speeds are about 1.5 m s −1 stronger than they are in the 3DPBL. As was 460 the case in the other two stabilities, MYNN still predicts stronger losses due to internal wakes. The unstable MYNN capacity factor losses (27.7 pp-29.6 pp) are substantially stronger than the unstable 3DPBL losses (10.4 pp-12.1 pp) because power output in Region II of the power curve is more sensitive at larger values of wind speed (due to the dependency on wind speed cubed). However, unlike in stable and neutral conditions, in unstable conditions MYNN has an overall higher average capacity factor, as its substantially stronger NWF winds offset its stronger wakes.

465
Paralleling the behavior of hub-height wind speed deficits, the inclusion of explicitly generated TKE consistently leads to smaller wake-induced power losses. The addition of TKE led to a reduction of hub-height wind speed deficits (Sec. 3.2).  As such, the losses in the 100TKE simulations are approximately 1 pp-2 pp smaller than they are in the 0TKE simulations, regardless of stability or PBL scheme.
In the end, these power calculations emphasize that power losses are a function of wind speeds coming into a plant, wake 470 recovery rate, and TKE addition. We stress that these idealized simulations have been carried out in pseudo-steady conditions, and they occur for one set of geostrophic winds in one part of the power curve. To better predict the cumulative non-linear interactions of the effects of these parameters on losses at a real-world location, we simulate a month-long case study in the U.S. mid-Atlantic coast. The interplay between incoming hub-height wind speeds and wake recovery rates becomes much muddier in a time-varying environment.

4 Results: Mid-Atlantic Case Study
In this section, we compare wind speed deficits in wakes produced with MYNN and the 3DPBL in August 2020 in the mid-Atlantic. As described in Sec. 2.4, we run three categories of simulations-NWF, Vineyard Wind 1 only, and all the lease areas.
We pay special attention to Vineyard Wind 1 because our set of simulations allows us to differentiate between internal and external waking at this site. Finding only minor TKE differences between MYNN and the 3DPBL in the idealized simulations, 480 we omit TKE analysis at this site and focus instead on wind speeds and power production.

Turbine-Free Winds
Before analyzing wake effects in the mid-Atlantic, we first examine wind speeds in turbine-free simulations. Specifically, we calculate average profiles at the Vineyard Wind 1 centroid. We classify each 5-minute interval output as stable, neutral, or unstable. While a number of metrics can be used to classify atmospheric stability (e.g., bulk Richardson number, Obukhov 485 length), we use WRF-predicted surface heat fluxes at the Vineyard Wind 1 centroid to facilitate comparison to the idealized simulations. We define stable conditions as having a heat flux less than -5 W m −2 , unstable conditions as having a heat flux greater than 5 W m −2 , and neutral conditions as in between. We designate these thresholds so they overlap with the stability metrics used in the idealized simulations. We emphasize that offshore heat fluxes in this domain are significantly weaker than typical heat fluxes observed on land. As such, while we refer to atmospheric states as "stable" and "unstable" in this offshore 490 environment, these states can resemble the atmosphere under near-neutral stratification. MYNN and the 3DPBL spend a similar percentage of the month under stable conditions (35% and 33%, respectively), whereas MYNN shows more frequent neutral conditions than the 3DPBL (50% versus 40%) and, conversely, MYNN shows less frequent unstable conditions than the 3DPBL (15% versus 27%).
The average NWF wind speed profiles at Vineyard Wind 1 share characteristics with wind speed profiles in the idealized 495 simulations (Fig. 9)  Vineyard Wind profiles are similar across PBL schemes, as was the case in the idealized simulations, but, notably, the 3DPBL now shows slightly faster hub-height wind speeds.
We also characterize the distribution of hub-height NWF wind directions and wind speeds using a wind rose at the Vineyard 505 Wind centroid (Fig. 10). Under stable and neutral conditions, winds are predominantly out of the southwest, whereas under unstable conditions, they are out of the northeast and presumably influenced by land, Martha's Vineyard. Thus, the wakes under unstable conditions will point to the southwest, opposite to those under neutral and stable conditions, which will be directed to the northeast. In all three stability conditions, the distribution of wind directions is fairly narrow, which facilitates the appearance of wakes in time-averaged visualizations. MYNN shows a greater prevalence of winds at rated power (Region III) under stable and unstable conditions. In neutral conditions, both MYNN and 3DPBL winds are predominantly in Region II; therefore, we expect less wind resource (larger wake effects on power production) in this stability condition, at least for the month of August.

Average Hub-Height Wind Speed Deficits
We calculate the average hub-height wind speed deficits by stability in the mid-Atlantic (Fig. 11). We time average the hub-515 height wind speeds in the NWF, VW-ONLY, and LEASE simulations, categorizing the stability of each 5-min period based off heat fluxes at the centroid of the NWF Vineyard Wind simulations and using the same timestamps for all three wind turbine cases for each PBL scheme. The distribution of wind direction is narrow in each stability and, as such, we do not additionally filter by wind direction in the wake analysis as might be done at a site with more variable wind directions.
Monthly averaged internal mid-Atlantic wakes, quantified at Vineyard Wind 1, are typically sensitive to both stability as 520 well as PBL scheme (Fig. 11a- internal waking than the 3DPBL (-1.02 m s −1 versus 1.06 m s −1 , or 18.5% versus 18.7%). We also note that even though average unstable hub-height NWF wind speeds are 2.0-4.3 m s −1 stronger than their neutral counterparts, the average absolute unstable MYNN internal waking is only 0.1-0.2 m s −1 stronger than under neutral conditions. Offshore wakes under unstable conditions in these real simulations erode much more quickly than they would when forced under neutral or stable stratification under the same winds, even when the surface heating is weak. This effect can also be seen by looking at relative waking 530 effects-relative neutral internal wakes are about 18.6%, whereas they are 12.3% in unstable conditions.
The different PBL schemes lead to different internal waking, although it is more difficult to characterize the exact mechanisms in the time-varying mid-Atlantic simulations than in the pseudo-steady-state idealized simulations. The differences in hub-height NWF wind speed play a key role. In the mid-Atlantic simulations, when MYNN has a stronger NWF hub-height wind speed, it produces a stronger averaged internal wake (as compared with stable and unstable conditions). The NWF wind 535 speed differences are smallest under neutral conditions in which the 3DPBL is only ∼0.25 m s −1 stronger than MYNN at hub-height. In these conditions, it is easier to observe the impact of different turbulent mixing formulations. In the idealized neutral simulations, MYNN had slightly stronger NWF wind speeds (coincidentally, ∼0.25 m s −1 stronger), and its wakes were also stronger (also by ∼0.25 m s −1 ). In the real simulations, 3DPBL NWF winds are slightly stronger, but the internal wakes have a near-equal magnitude between the PBL schemes (differences <0.05 m s −1 ). Thus, the 3DPBL's faster internal 540 recovery rate observed in the idealized neutral simulations appears to persist through to the mid-Atlantic simulations as well.
Aside from differences in NWF wind speeds and turbulent mixing approaches, the two PBL schemes also may lead to different wake magnitudes due to effects from time-varying winds as well as interactions with other physics parameterizations that were not activated in the idealized simulations.
Monthly averaged external wake propagation is particularly sensitive to stability and only weakly sensitive to the PBL 545 scheme choice (Fig. 11d-f,j-l). We characterize external wake propagation in the mid-Atlantic by highlighting the 0.5 m s −1 contour and neglecting the e-folding distance, because the 0.5 m s −1 contour is a sufficiently useful metric when all internal wakes are substantially stronger than 0.5 m s −1 . As expected, external wakes propagate furthest under stable conditions. In the stable VW-ONLY cases, they propagate 51 km east of the easternmost point in Vineyard Wind with MYNN and 38 km east with the 3DPBL (Fig. 11a, g). In the LEASE cases, these external wakes grow in size, although a characteristic length is more 550 difficult to quantify due to their irregular shape. In neutral conditions, the 0.5 m s −1 deficit remains close to the plant with both PBL schemes because of the weaker overall wind speeds. While the external wakes under unstable condition also remain fairly close to the plants (again, due to overall weaker wind speeds), they do propagate slightly further in the MYNN simulations (due to overall stronger wind speeds).
Additionally, we quantify the impact of external wakes on hub-height wind speed deficits by focusing on wind speed deficits 555 inside Vineyard Wind 1. We calculate the monthly averaged external wake effect as the average internal wake magnitude in the VW-ONLY simulations subtracted from the perceived internal wake at Vineyard Wind in the LEASE simulations (e.g., Fig. 11a subtracted from Fig. 11d). Unsurprisingly, the largest external wake effects are observed under stable conditions. MYNN has a stronger external wake effect (0.78 m s −1 or 6.4 pp) than the 3DPBL (0.56 m s −1 or 5 pp) in these conditions. Interestingly, under unstable conditions, the 3DPBL has a stronger wake effect even though its internal waking is weaker than MYNN's.

Impact on Power Production
We quantify power production at Vineyard Wind 1 in the VW-ONLY simulations for each of the stability conditions ( Fig. 12a- MYNN's stronger NWF winds. This large discrepancy stems from the differences in wind speed distributions (Fig. 10c, f).
Depending on the PBL scheme, power losses associated with internal wakes disagree on the order of 1-2 percentage points. 575 We calculate the expected unwaked power production in the idealized simulations by convolving the time-varying hub-height wind speeds with the power curve, and we calculate internal wake power loss with reference to this value. Under stable conditions, MYNN predicts slightly stronger losses than the 3DPBL (26.2 pp loss versus 24.9 pp loss). Under unstable conditions, we see the same general behavior-average MYNN losses are 12.4 pp and average 3DPBL losses are 10.3 pp. Even though winds were much weaker overall under neutral conditions, we also see an approximate 1-pp difference in power loss. This stronger internal waking across all stability conditions. The mid-Atlantic simulations also highlight that even though average capacity factors incorporating wake effects may be drastically different between PBL schemes (e.g., 56.1% for MYNN versus 38.2% for the 3DPBL), the internal wake losses can still be similar.
External wakes can amplify losses ( Fig. 12d- factors that are 1.6 pp-2.8 pp less than in MYNN-or about 4.7%-7.8% less total power generated. We additionally visualize timeseries of wind speeds and power production to underscore the key role that NWF wind speed plays on power production (Fig. 13). We highlight NWF hub-height winds, capacity factor, and wake effects for one week of the simulation. In general, when NWF hub-height winds at the Vineyard Wind 1 centroid are stronger for a given PBL scheme ( Fig. 13a), that PBL scheme also produces more instantaneous power (Fig.13c). This amplification is particularly true when 600 NWF wind speeds are in the vicinity of the rated wind speed (between Region II and Region III), as occurs on August 23. This pattern also persists when winds are weaker (as they are between August 24-25), although it is less prominent because the wakes are also weaker. Differences in wind speed are less important when NWF winds exceed the rated speed (as on August 26), as internal wakes in the VW-ONLY simulation do not reduce power output. Even then, external wakes can further weaken winds inside Vineyard Wind and reduce power output, as they momentarily do on August 26. Regarding external wakes, even 605 minor differences in NWF wind direction between PBL schemes (< 5 • , Fig. 13b) can produce instantaneous external wake effects that differ by a dozen pp or greater (Fig. 13e), as occurs on August 23.
While faster wind speeds tend to lead to a larger capacity factor for a given PBL, faster wind speeds do not necessarily cause stronger internal or external waking. During the first red-highlighted period of Fig. 13, MYNN has faster NWF wind speeds but weaker internal waking. This disparity could be explained by subtle differences in NWF wind direction, significant 610 spatial variability of wind directions within Vineyard Wind 1 (rendering the single-point measurement of wind direction less informative), or differences in mixing for each PBL scheme, among other factors.
In the end, these timeseries figures show that the two PBL schemes can predict capacity factors that differ by dozen of pp in a given moment. Thus, when seeking to characterize power production uncertainty, it may be even more beneficial to vary PBL schemes for short-term wind forecasts than for long-term wind resource assessments.   In the end, either PBL scheme could lead to stronger total power production in the idealized plants, depending on the relative effects of NWF winds and wake recovery rates. In general, however, MYNN produced faster winds and stronger wakes. In the stable and neutral 100TKE simulations, the 3DPBL predicted average capacity factors that were about 1-2 percentage points higher than with MYNN (or produced about 3%-6% more overall power). However, due to the substantially faster NWF 640 winds under unstable conditions, MYNN predicted an average capacity factor that is about 13 percentage points higher than with the 3DPBL (or produced about 35% more overall power). While either PBL scheme could produce a higher capacity factor depending on stability, MYNN consistently had stronger losses due to internal wakes. MYNN capacity factor losses in the 100TKE simulations were 9.9-27.7 percentage points; they were 6.6-10.4 percentage points for the 3DPBL. When explicit TKE addition was turned off, power losses due to internal wakes consistently increases by about 1-2 percentage 645 points, independent of stability and PBL scheme.
Wind speed deficits in the real-world mid-Atlantic simulations also depended on the relative importance of NWF winds and wake recovery rates, although disentangling their relative importance in these time-varying fields became more difficult. While the two PBL schemes produced substantially different capacity factors, they predicted relatively similar power losses due to internal wakes, differing by only 1-2 percentage points. Average internal losses at Vineyard Wind 1 were about 25.5 665 ± 0.6 percentage points under stable conditions, 10 ± 0.5 percentage points under neutral conditions, and 11 ± 1 percentage point under unstable conditions. Average external wakes showed slightly more variability between PBL schemes but, in the end, they further reduced average capacity factors at Vineyard Wind 1 by about 6-8 percentage points under stable conditions, 2 percentage points under neutral conditions, and 3-4 percentage points under unstable conditions. While the two PBL schemes lead to different monthly averaged assessments of power production and wake losses, these differences could be substantially 670 larger at a given moment. At times, instantaneous capacity factors differed by more than 20 percentage points between the two PBL schemes, and internal wake effects, as well as external wake effects, also showed momentary discrepancies of a similar magnitude. Thus, short-term forecasts of wind power production may be especially sensitive to PBL scheme.
In the end, this analysis suggests wind farm parameterizations are indeed sensitive to the choice of PBL scheme. In general (but not always), MYNN predicts stronger wind speed deficits inside of a plant than the 3DPBL, and MYNN tends to predict 675 stronger power losses due to internal wakes. But, MYNN also tends to predict stronger plant inflow wind speeds. Because of these competing effects (inflow wind speed vs wake dissipation), a generalization is not possible for which scheme consistently predicts larger or smaller power outputs.
Due to this sensitivity, we recommend that future wind energy planning studies that examine mesoscale model sensitivity consider varying the PBL scheme, along with other model inputs that have been established in literature, such as grid resolution, 680 magnitude of explicit TKE addition, and the choice of wind farm parameterization (Fischereit et al., 2021). While it is difficult to judge the sensitivity of PBL scheme choice relative to most of the other inputs, we note that changing the PBL scheme had a larger impact than turning explicit TKE generation on or off in the idealized simulations. While this analysis centered on Further, because the 3DPBL scheme is intended for mesoscale simulations approaching the "gray zone" where horizontal grid resolution is100-1000m, these 2-km resolution simulations are at the edge of where its representation of horizontal turbulent eddies are necessary. While few simulations with the wind farm parameterization have delved into this grey zone, it would be interesting to extend these simulations to finer resolution within the grey zone to determine optimal choices of scheme. Such extensions should compare simulated wakes to observed wakes, the next step in this work.

690
Code and data availability. Namelists for all simulations, time-averaged idealized WRF data, WRF code modifications, and analysis notebooks to reproduce all figures can be found on Zenodo (https://doi.org/10.5281/zenodo.5565399). For convenience, much of the same material has also been uploaded to GitHub (https://github.com/rybchuk/wfp-3dpbl-sensitivity). Data for the month-long 3DPBL NWF and LEASE mid-Atlantic simulations can be found by looking inside of year-long simulation data stored on the AWS Registry of Open Data at a link that will be updated during the peer review process. Data for the month-long 3DPBL VW-ONLY simulation as well the MYNN 695 simulations is stored on the University of Colorado's PetaLibrary and is available upon request.
Author contributions.