We utilize an eddy-resolving oceanic physical-biogeochemical model (ROMS-BEC) of the California Current System to generate average daily physical and biogeochemical conditions from 2000–2005 (Deutsch et al., 2021; Renault et al., 2021). ROMS simulates the ocean currents by solving the hydrostatic, primitive equations (Shchepetkin and McWilliams, 2005) in a terrain-following coordinate system and uses a K-Profile Parameterization (KPP) for vertical mixing (Large et al., 1994). These ocean currents partially control the coupled biogeochemical variability simulated online with BEC (described in Deutsch et al., 2021). BEC simulates distinct nutrient cycles (nitrogen, silicic acid, phosphate, and iron), organic matter, dissolved oxygen, and carbonate chemistry.

The ROMS-BEC simulations comprise a family of one-way nested solutions that successively increase resolution and decrease spatial extent (Mason et al., 2010). For the present analysis, the horizontal resolution is 1 km with 60 vertical levels, nested down from a 4-km parent solution (see Kessouri et al., 2020; Renault et al., 2021). This solution contains tidal forcing and synoptically variable atmospheric forcing prescribed by a 6-km resolution Weather Research and Forecast Model (WRF, Skamarock et al., 2008). Relevant ROMS-BEC parameters utilized by the macroalgal growth model include nitrate, ammonium, dissolved organic nitrogen, temperature, photosynthetically active radiation (PAR), concentration of chlorophyll-a (for light attenuation due to phytoplankton shading), and seawater velocities. Wave characteristics for the corresponding historical time period were downloaded from observational buoy data available from the National Data Buoy Center (NDBC; Station 46053).

The ROMS-BEC inputs to the macroalgal growth model exhibit both seasonal and higher-frequency variability (Figure 2), indicative of the complex flows and inherent nonlinearity of the coupled biogeochemistry. The ROMS-BEC simulation has a horizontal resolution of 1 km that fully resolves mesoscale variability and partially captures submesoscale variability (McWilliams, 2016); both regimes will impact on biogeochemistry (Kessouri et al., 2020). Specific to the shelf, higher resolution (e.g., on the order of 100 m) simulations produce small-scale flow patterns in the nearshore that can impact material fluxes and thus biogeochemistry (Dauhajre et al., 2019), however, computational costs (time and storage) limit production of multi-year, higher resolution simulations. Despite these computational constraints, we expect the 1-km simulation to adequately capture the seasonality of the inputs to the macroalgal growth model.

High-resolution Large-Eddy Simulations (LES) are used to study flow within and around the macroalgal farm, capturing motions in the range of scales between one and 100 m. In LES, most of the unsteady three-dimensional turbulent eddies are explicitly resolved on the numerical grid, and only the small-scale eddies that are more universal need to be parameterized (Pope, 2000). Because it requires fewer assumptions about the flow, this approach is particularly well suited for studying complex flows in which the main features are not well understood. LES has been used in studies of ocean turbulence since the mid 1990's (Skyllingstad and Denbo, 1995; McWilliams et al., 1997). We adopt the usual approach of solving the wave-averaged Navier-Stokes equations, in which surface wave motions are not explicitly resolved but their impact in the flow is represented via Stokes drift, together with the Boussinesq approximation. The LES approach and model setup employed here are described in detail in Yan et al. (2021), and only a brief description is provided below. For more details on the application of LES to ocean flow, please refer to the review paper by Chamecki et al. (2019).

The set of equations used in this study consists of the continuity equation (conservation of mass) under the Boussinesq approximation (1)∇·u=0, the wave-averaged Navier-Stokes equations (2)∂u∂t+u·∇u=-∇Π+fez×ug︸I-fez×u︸II+us×(ζ+fez)︸III                +(1-ρ/ρ0)gez︸IV-FD︸V+∇·τ︸VI, and the advection-diffusion equation for seawater density fluctuations (caused by temperature fluctuations) (3)∂ρ∂t+(u+us)·∇ρ=∇·τρ. In this system of equations, seawater density ρ, Eulerian velocity u = (u, v, w), and generalized pressure Π are evolved on the three-dimensional numerical grid represented using cartesian coordinates x = (x, y, z), with z being the vertical direction and e_z the vertical unit vector. In addition, ρ₀ is the reference density, f is the Coriolis frequency, g is the gravitational acceleration, and u_s and u_g represent the Stokes drift velocity and the mean current in geostrophic balance, respectively. The effect of the small scales (i.e., the subgrid scales; SGS) on the resolved scales is represented by the SGS momentum flux tensor τ and the buoyancy flux vector τ_ρ. Molecular viscosity is neglected on the basis of the high-Reynolds number flows considered here.

The terms on the right-hand side of Equation (2) are the pressure-gradient force (I), the Coriolis force (II), the forces resulting from the wave-averaging procedure (III), the buoyant force from the Boussinesq approximation (IV), the drag force exerted by the macroalgal farm (V), and the SGS force (VI). The pressure-gradient force is split into a fluctuating component (−∇Π) and a mean component that is used to drive the mean flow (fe_z × u_g, where u_g is an imposed geostrophic velocity that in practice determines the strength and direction of the mean pressure-gradient force). The first of the two terms including the Stokes drift is the vortex force us×ζ~ with ζ~=∇×u~ being the vorticity vector (Craik and Leibovich, 1976), which is key to driving Langmuir turbulence in the upper ocean (Leibovich, 1983), and the second term represents the interaction between Stokes drift and Coriolis. For simplicity, we only consider steady monochromatic waves, so that the Stokes drift velocity is given by (4)us(z)=Usexp(2kz)ew. In Equation (4), Us=ωkaw2 is the Stokes drift velocity at the surface, k is the wavenumber, e_w is a unit vector in the direction of the wave propagation, and a_w and ω are the wave amplitude and angular frequency.

Solution of Equations (1)–(3) require closure models for the SGS terms and the drag force exerted by the farm. The SGS momentum flux is parameterized using the scale-dependent dynamic Smagorinsky model (Bou-Zeid et al., 2005), and the SGS density flux is parameterized using a gradient-diffusion model with diffusivity obtained from the SGS viscosity and a constant Schimdt number Sc = 0.4 (Yang et al., 2015). Finally, as the numerical grid is too coarse to resolve the geometry of the macroalgae fronds and stipes, we represent their effect on the flow through a body force that represents the total drag within a grid cell given by, (5)FD=12CDaP·(|u~|u~). Here, C_D = 0.0148 is the drag coefficient according to the field experiments by Utter and Denny (1996), a is the total frond surface area per unit volume of space (obtained by conversion of the algal biomass; Supplementary Table 1), and the projection coefficient tensor P acts to yield the effective surface area facing each direction (in the present case we use P = (1/2)I, where I is the 3 × 3 identity matrix). The motion of macroalgae elements in response to the flow can be estimated based on their geometry, modulus of elasticity and buoyancy (Luhar and Nepf, 2011; Henderson, 2019). As a first approximation, we assume that the macroalgae move with the wave orbital wave velocity and thus pose no impact on the wave field (Yan et al., 2021).

Thus, in the modeling approach adopted here, the ocean state experienced by the farm is mostly determined by the surface fluxes (wind stress u_∗ and surface buoyancy flux B₀), the mean current (u_g), the Stokes drift profile (imposed by specifying wave parameters k, a_w, and ω), and the stratification in the pycnocline (represented by a constant temperature gradient). These forcings evolve continuously on diurnal and seasonal time scales, producing a wide range of ocean conditions. LES is not suitable to produce a continuous year-long simulation of evolving ocean conditions due to computational limitations. In the present work, we adopt a simplified approach in which we use LES to study a few “typical” conditions with steady forcing (i.e., no temporal changes in the ocean state) and then use ROMS-BEC results to modulate the flow features uncovered in these few LES runs (see below for more details).

We consider two cases with different surface forcing leading to distinct upper mixed-layer depths: (1) a deep, upper-boundary-layer case characteristic of winter and (2) a shallower, upper-boundary-layer case typical of summer. In both cases, we simulate neutral boundary layers with no incoming or outgoing buoyancy flux at the surface (i.e., B₀ = 0). Fluid below the mixed layer is stably stratified with a constant initial potential temperature gradient dθ/dz = 0.01 °C m⁻¹. In the deep boundary layer case, the wind stress at the surface is u*=6.1×10-3 m s-1, corresponding to a wind speed at 10-m height above the surface of U10=5 m s-1. The surface wave has a wavelength of λ = 60 m (k = 2π/λ) and a wave amplitude of a_w = 0.8 m. The resulting Stokes drift at the surface is Us=0.068 m s-1, yielding a turbulent Langmuir number Lat=(u*/Us)1/2=0.3 that is representative of equilibrium wind-wave conditions in the open ocean (Belcher et al., 2012). This results in a boundary-layer depth z_i ≈ 25 m. For the shallow, upper-boundary-layer case, we reduce the surface forcing to u*=3.05×10-3 m s-1 and Us=0.034 m s-1, resulting in a shallower upper boundary layer with z_i ≈ 15 m. In both cases, a steady current ug=0.2 m s-1 is also imposed to represent the effect of ocean currents or mesoscale eddies (assuming that mesoscale spatial and temporal variability is negligible in the LES domain). This current speed occurs through the year, but is more characteristic of the summer (Figure 2E). For the sake of simplicity, we assume that the wind, the waves, and the mean current are all aligned in the same direction. More details about the numerical model and the simulation setup can be found in Yan et al. (2021), and the geometry of the farm used in the LES runs is described in Section 2.4.

The macroalgal growth model (MAG) is developed for Macrocystis pyrifera. Existing macroalgal growth models informed the conceptual framework of the model presented here (Jackson, 1987; Solidoro et al., 1997; Broch and Slagstad, 2012; Hadley et al., 2015). We have compared our model selection to alternative forms of the growth model and found that our calculations of biomass are robust to model selection (Appendix A). Nitrogen is tracked in six pools between seawater and macroalgae: nitrate (NO₃), ammonium (NH₄), dissolved organic nitrogen (DON), particulate organic nitrogen (PON), stored nitrogen in macroalgae (N_s), and fixed nitrogen in macroalgae (N_f) (Figure 3A). The currency of nitrogen has been selected as this is the nutrient limiting productivity for M. pyrifera along the U.S. West Coast (Zimmerman and Kremer, 1984, 1986). For conversion to carbon units, M. pyrifera has an average C:N ratio of 16.6, but this value varies greatly in natural kelp, between 6 and 48, given macroalgae's ability to store intracellular pools of both carbon and nitrogen (Rassweiler et al., 2018). Partitioning between pools of internal nitrogen (N_s and N_f) is a two-step uptake process where nitrogen is first stored in intracellular pools and then assimilated into the macroalgae's cellular structure and is adapted from Hadley et al. (2015). There is a minimum amount of nitrogen required, Q_min (mg N g-dry^-1), and the nutrient quota, Q, defines the relative nitrogen status of macroalgae. (6)Q=Qmin(1+NsNf) Novel aspects of the MAG model developed here are the accounting of individual fronds, the inclusion of demographic rates of fronds such as initiation, age, growth, and senescence, and the upward growth of macroalgae through the water column. Holdfasts and sporophylls are not included in the model as they contribute little to uptake and growth (Gerard, 1982; Arnold and Manley, 1985). Fronds are categorized as subsurface, water-column, and canopy (following Rassweiler et al., 2018, depicted in Figure 1). “Subsurface” fronds are those that have not reached the surface. Fronds that have reached the surface are treated in two sections: the “water-column” section is the portion of frond that is underwater stretching from the cultivation depth to the sea surface and the “canopy” section is the portion of the frond that floats at the sea surface. Conversion between frond nitrogen and other allometric parameters is based on relationships that have been established for M. pyrifera (Supplementary Table 1). The conversion between nitrogen content and biomass is calculated as: (7)B=NfQmin Biological parameters used by MAG are provided in Table 1. New fronds are initiated, S_i, as a function of the average nutrient quota, Q, of existing fronds given that light at depth is greater than the compensating light irradiance. New fronds are initiated with N_s,i and N_f,i that is equivalent to 1-m frond with an average Q of existing fronds. Frond senescence is age dependent with a maximum life span of Age_max (Rodriguez et al., 2013).

Macroalgal growth is modeled as a two-step process. Nitrogen is taken up from seawater into an internal reserve of stored nitrogen (N_s; mg N m^-3). N_s is converted into fixed nitrogen via growth (N_f; mg N m^-3). The change in N_s over time is a result of uptake of nitrogen from seawater (U; mg N g-dry^-1 d^-1), growth (G; d^-1), exudation (E; d^-1), mortality (M; d^-1), harvest (H), and frond initiation [S_i; Equation (8)]. (8)∂Ns∂t=U×B-(G+E+M+H)×Ns+Si×Ns,i Since M. pyrifera is buoyant, fronds extend upward. N_s is vertically resolved along the length of the frond at a resolution of dz = 1 m. Translocation of N_s is a physiological phenomenon in which N_s is redistributed through a network of sieve tubes on the scale of multiple meters per day (Parker, 1965; Schmitz and Lobban, 1976). MAG accounts for this process by redistributing N_s along the length of the frond relative to the vertical distribution of fixed nitrogen. This maintains a constant Q along the length of the frond. Exudation is included here as the chemical composition of exudates has been shown to include nitrogen as well as carbon (Chen et al., 2020).

The change in N_f over time is due to growth (G; d^-1), mortality (M; d^-1), harvest (H), and frond initiation [S_i; Equation (9)]. (9)∂Nf∂t=G×Ns-(M+H)×Nf+Si×Nf,i The entire frond contributes to growth, but new N_f is added at the tip of the frond. N_f moves into new depth bins when a capacity term (N_f-to-length; mg N m^-1 frond) in the existing depth bin is reached. This capacity term is defined as the maximum amount of fixed nitrogen that can be contained per meter of frond length (Supplementary Table 1). In the canopy, N_f growth ceases as total frond height approaches maximum frond height (h_max). N_f in the canopy is redistributed horizontally based on the average length of the canopy portion of the fronds using a uniformly weighted smoothing function.

Nitrogen sources for M. pyrifera are ammonium, nitrate and urea (Haines and Wheeler, 1978; Gerard, 1982; Harrison and Hurd, 2001; Smith et al., 2018). These nutrients are transported in seawater with relevant sink and source terms [Equations (10)–(12)]. The concentration of urea is calculated as 20% of DON (Sipler and Bronk, 2015). (10)∂NO3∂t=Diffusion-Advection-UNO3×B (11)∂NH4∂t=Diffusion-Advection-UNH4×B              +τDON×DON (12)∂DON∂t=Diffusion-Advection-Uurea×B            +(E+M)×Ns-τDON×DON where, U is the nutrient-specific uptake rate, and τ_DON is the remineralization rate of DON back to NH₄. Details of 3-D nutrient transport [Equation (29)], including advection and diffusions terms, are provided in Section 2.4. Macroalgae also contribute to a particulate organic nitrogen (PON) pool via mortality: (13)∂PON∂t=M×Nf The total uptake rate, U, is the sum of the uptake rates for nitrate, ammonium, and urea (Haines and Wheeler, 1978), which is kinetic and mass-transfer limited (Stevens and Hurd, 1997), as well as dependent on the nutrient quota, Q (Kopczak, 1994). Nitrogen transport occurs across the membranes of cells at the blade with negligible contribution from stipe (Gerard, 1982). Uptake rates are computed as surface area fluxes and converted to dry weight (DW) fluxes. (14)U=∑i=13Umax,i[f(Ci,|u|,Tw)×f(Q)]×AreaDW×2BbladeBfrond where, C_i is the concentration of the constituent nutrient in seawater (i(1:3) = [nitrate, ammonium, urea]), |u| and T_w are the velocity magnitude and wave period, respectively. The fraction of the frond biomass (B_frond) that is blade (B_blade) is based on an empirical relationship of blade biomass with fractional frond height (Nyman et al., 1993). The limitation of uptake due to kinetic and mass-transfer limitations is modeled similar to Stevens and Hurd (1997) as: (15)f(Ci,|u|,Tw)=CiKm,i(CiKm,i+12(γ+γ2+4CiKm,i)) where (16)γ=1+Vmax,iβKm,i-CiKm,i

K_m,i is the half saturation constant and V_max,i is the maximum uptake rate for nutrient i (Table 1). These parameters have been estimated for each nutrient from both laboratory and in situ time-course experiments. β in the above formulation is: (17)β=DδD+4δDTw∑n=1∞(1-exp(-Dn2π2Tw2δD2)n2π2) The first term represents mass-transfer limitation due to the presence of a diffusion boundary layer, δ_D, and the second term represents the stripping away of the diffusion boundary layer twice each T_w. D is the molecular diffusion coefficient (Yuan-Hui and Gregory, 1974). δ_D is modeled as an empirical function of the magnitude velocity with kinematic viscosity, v (Stevens and Hurd, 1997). (18)δD=10v0.33|u| The limitation of uptake due to internal nutrient reserves is based on the observation that macroalgae have nitrogen storage capacities that reflect past nutrient supplies (Hurd et al., 2014). Luxury uptake is formulated so that uptake decreases linearly with increasing Q as: (19)f(Q)=Qmax-QQmax-Qmin The maximum growth rate, μ_max is limited by temperature (T), availability of photosynthetically active radiation (PAR), and internal nutrient reserves (Q). (20)G=μmax[g(T)×g(PAR)×g(Q)]

M. pyrifera sporophytes exist across a broad range of temperature, but do not survive well in very cold or very warm waters (Schiel and Foster, 2015). In the field, growth rates have been observed to be independent of temperature below 21 °C (Zimmerman and Kremer, 1986). The temperature limitation on growth is expressed as a piecewise formulation (Table 1): (21)g(T)={T/Tmin,T<Tmin1,Tmin≤T<Tmax-1/3×T+7.7,Tmax≤T<Tlim0,T≥Tlim} The effect of PAR on growth is derived from field-based transplant experiments of juvenile M. pyrifera (Dean and Jacobsen, 1984). (22)g(PAR)={0,PAR<PARc1-exp[kPAR(PAR-PARc)],PAR≥PARc}

k_PAR is the light-limited growth constant (Table 1). PAR_c is the value at which growth ceases, known as the light-limited compensating irradiance (Table 1). PAR is depth-resolved and attenuates due to optical properties of seawater (0.0384 m^-1), shading by phytoplankton (0.0138 m² mg^-1 chl-a), and self-shading by macroalgae (Table 1, Lorenzen, 1972).

Growth rate is dependent on nutritional history, increasing linearly with increasing Q as: (23)g(Q)= Q-QminQmax-Qmin There are both biological and physical contributors to mortality (Table 1). d_wave is wave-dependent loss and is a function of the significant wave height, H_s (Rodriguez et al., 2013). d_blade is continual loss of blades due to deterioration (Rassweiler et al., 2018), applied to the fraction of frond that is blade. d_frond is the fractional rate of loss once a frond has reached Age_max. We do not include all forms of potential mortality, such as grazing and pests, as it is unclear how this source of mortality will translate to suspended farm settings. (24)M={dwave+dblade×BbladeBfrond,Age<Agemaxdfrond,Age≥Agemax} Harvesting of Marocysitis pyrifera entails removing the canopy portion of fronds while the subsurface fronds stay intact. The subsurface portion of cut canopy fronds quickly senesce at a rate of d_frond, while uncut fronds continue to grow, regenerating a new canopy.

The simulated farm is located 3 km from the coast, 6 km from the nearest harbor, and directly south of the Mohawk Reef Kelp Forest in Santa Barbara, CA (119.733 W 34.368 N; Figure 2A). The design is 400 x 400 m (16 hectares) with longlines extending the length of the farm (x-direction), repeated every 26 m (y-direction), and oriented parallel to the dominant flow direction. Growth lines, where kelp is seeded, are attached perpendicular to the longline every other meter and extend 8 m in length. Kelp is cultivated at 20 m below the sea surface (z_cult) in 60 m water depth, such that longlines are suspended in the water column via an anchor-and-tether design (Figure 1B). A set of LES simulations are first performed to generate modified flow features that are then represented in the full-farm simulations. The full-farm simulations integrate ROMS-BEC upstream conditions, LES within-farm modified flows, and the MAG model to predict spatially resolved kelp biomass.

For the LES runs, the longlines are assumed to align with the mean current and the wind/wave propagation direction. Even though the angle between the farm longlines and the flow is expected to have a large impact in the flow field, this investigation is outside the scope of the present paper. Note that with the chosen cultivation depth, the farm is completely within the upper, well-mixed boundary layer in the deep layer simulation (winter conditions) but it extends below the upper mixed layer in the shallow case (summer conditions). For each of the two ocean conditions we perform three LES simulations representing different growth stages of the canopy (total of six simulations): (i) a no-drag simulation representing initial stages of growth in which the farm poses very little drag, (ii) a subsurface stage that represents an intermediate growth stage with kelp extending 17 m above the farm base, and (iii) a fully grown stage with a large canopy on the surface (Supplementary Table 3). Selection of which LES simulation to use during the year-long full-farm simulations is based on the upstream ROMS-BEC depth of the upper mixed layer. If the upper mixed layer depth is less than the cultivation depth, results from the shallow simulation are used, and if the upper mixed layer depth is greater than the cultivation depth, results from the deep simulation are used. To select which LES results to use based on the growth stage of the farm, we evaluate the canopy height across the farm. Results from the no-drag LES simulation are selected when the farm-averaged maximum frond height is less than 10 m, results from the subsurface LES simulation are selected when this value is between 10 and 24 m, and results from the fully-grown LES simulation are selected when this value is greater than 24 m.

As it will be shown in Section 3.1, LES results reveal two main flow features produced by the presence of the farm that must be represented in the MAG model: secondary flow structures in the time-averaged flow (termed attached Langmuir circulations) and modulations of the turbulence intensity. While the former are represented by modification of the ROMS advection velocities, the latter are included as modulations of the KPP eddy diffusivity. Note that the ROMS provides a time evolving profile of the horizontal velocity field u_ROMS(z, t) upstream of the farm.

To represent the attached Langmuir circulations in the full-farm simulations, we define a 3-D time-averaged velocity field from the LES, u¯(x,y,z), and an associated velocity scale based on the mean velocity component parallel to the longlines upstream of the farm (i.e., at x = 0) given by U0=〈u¯〉y,z, where angle brackets indicate spatial averaging in the direction indicated and the average in z is performed within the entire vertical extent of the boundary layer. From this, a 3-D dimensionless velocity deficit is defined (25)ΔuLES(x,y,z)=u¯(x,y,z)-U0U0. This velocity deficit is then applied to the vertically-averaged (from the surface to the mixed layer depth) time-varying ROMS flow U_ROMS(t) = 〈u_ROMS〉_z(t) to determine the velocity used in MAG (26)uMAG(x,y,z,t)=〈uROMS〉z(t)+UROMS(t)ΔuLES(x,y,z). In this way, the time evolution of the advection velocity experienced by the farm is determined by ROMS, but the spatial patterns within and around the farm (including the attached Langmuir cells) are represented by the velocity deficit obtained from the LES field.

To represent the effects of the farm on turbulence (and consequent enhanced nutrient transport), we estimate a vertical diffusivity from the LES runs using a modified version of the classic k-ϵ closure, where k represents the turbulent kinetic energy (TKE) and ϵ is its rate of dissipation (Pope, 2000). This type of closure is frequently used to represent mixing within terrestrial canopies (Katul et al., 2004). Because the focus is in modeling the vertical eddy diffusivity D_v, we use the vertical velocity variance σw2(x,y,z) instead of the full TKE as, (27)Dv,LES(x,y,z)=Cμσw4ϵ, where ϵ(x, y, z) is the 3-D field of the TKE dissipation rate estimated from the SGS dissipation in the LES. Here, C_μ = 0.2 is an empirical constant calibrated so that the model, when applied to the flow upstream from the farm, is in good agreement with the eddy diffusivity for Langmuir turbulence diagnosed by McWilliams et al. (1997). Similar to the strategy for the advection velocity, the diffusivity from LES is used to modulate a time evolving diffusivity from ROMS (28)Dv,MAG(x,y,z,t)=Dv,LES(x,y,z)〈Dv,ROMS(z,t)〉z〈Dv,LES(x=0,y,z)〉y,z.

For the full-farm simulations, four 3-D farm simulations are performed (Simulation A-D; Table 2). Each simulation differs in the implementation of the within-farm nutrient transport equation in MAG to identify which processes in the physical and nutrient environment impact farm performance [Equation (29)]. Simulations A-C do not include harvest of macroalgae, Simulation D does. (29)∂C∂t=-∂(uC)∂x︸I-∂(vC)∂y︸II-∂(wC)∂z︸III+∂∂x(Dh∂C∂x)︸IV            +∂∂y(Dh∂C∂y)︸V+∂∂z(Dv∂C∂z)︸VI±S︸VII where, C is the concentration of the nutrient being transported, u, v, and w are the velocities in the x, y, and z direction, D_h and D_v are horizontal and vertical diffusivity, and S represents the source and sink pathways specific to the nutrient being transported. D_h is 0.1 m² s^-1.

The transport equation for Simulation A includes horizontal transport terms (Equation (29) terms I, II, IV, and V) but vertical transport and sink terms are set to zero (Equation (29) terms III, VI, and VII). Horizontal advection is informed by the upstream depth-varying and time-varying ROMS-BEC conditions. Transport is not modified by LES results. This simulation is most similar to current macroalgal growth models that do not incorporate within-farm hydrodynamics. The transport equation for Simulation B includes horizontal transport terms and the source and sink term specific to the nutrient being transported (Equation (29) terms I, II, IV, V, and VII), but vertical transport terms are set to zero (Equation (29) terms III and VI). The horizontal advection terms are informed by upstream depth-varying and time-varying ROMS-BEC conditions. Simulation C and D implement the full transport equation (Equation (29) terms I–VII), while Simulation C does not include kelp harvest and Simulation D does. The advection and diffusion terms are informed by upstream ROMS-BEC boundary conditions and modified by LES results using Equations (26) and (28).

Full-farm simulations are initiated by seeding growth lines with biomass equivalent to a single 1-m frond at the beginning of the year and duration is one year. Spatial resolution is 1 m³. We also perform single-day full-farm simulations to demonstrate within-farm modified nutrient fields. Within-farm transport is solved at a resolution of 30 s. Determination of total kelp nutrient uptake is conditional; re-calculated when the maximum delta-N anywhere within the farm exceeds 20%, which varies from every 10 minutes to once per day depending on upstream conditions and stage of kelp. Kelp growth and mortality are solved once per day. For Simulation D, harvest is included and reported as metric tons of fresh weight per hectare (WMT ha^-1). The criteria to harvest is that there needs to be more then 2.5 WMT ha^-1 available in the canopy and biomass loss exceeds growth.

Figures 4, 5 show a summary of the results from the large eddy simulation for the full canopy farm in a deep upper-mixed layer (winter conditions). The drag forces exerted by the longlines of macroalgae on the flow (Figure 4A) cause a large reduction in the horizontal flow (Figure 4B), producing mean downwelling at the leading edge of the farm. The downwelling region extends well into the farm, giving place to upwelling near the trailing edge (Figure 4E). The spatial distribution of the drag caused by the longline spacing produces a pattern of alternating fast and slow currents (note the large negative deficit within the canopy rows in Figure 4B). This crosswise shear in the streamwise velocity (approximately parallel to the longlines) is associated with a stationary pattern of vertical vorticity which, in the presence of Stokes drift promoted by surface waves, leads to pairs of counter-rotating streamwise vortices within the farm (Yan et al., 2021). These vortices are seen by the alternating pattern of upwelling and downwelling shown in the mean vertical velocity [Figure 4D, and also in Yan et al. (2021) figures 6 and 12]. Because the dynamics behind the formation of these vortices is similar to the Craik-Leibovich type II instability (Leibovich, 1983), they were termed attached Langmuir cells in Yan et al. (2021). Note that the upwelling regions tend to occur within the macroalgal rows while the downwelling is mostly located between rows. In the present farm configuration, the attached Langmuir cells extend from the top to the bottom of the farm and are present in all four LES farm simulations used here (i.e., they occur for subsurface and full canopy cases both in shallow and deep mixed layers; see Supplementary Figures 1, 3, 5).

In addition to these changes in the mean flow, the drag forces exerted by the macroalgae impact turbulence within the farm. Here, we quantify this effect by calculating an eddy diffusivity based on the modified K-ϵ model (Equation 27). The vertical distribution of D_v upstream of the leading edge (Figures 5B,C) is consistent with that observed in Langmuir turbulence (McWilliams et al., 1997), with a peak in the upper half of the boundary layer. As the macroalgae plants reduce the turbulence levels and increase the energy dissipation within the canopy rows, the vertical diffusivity is found to be greater in the spacing between farm rows than that in the regions occupied by canopy (Figure 5). Contrary to the mean flow patterns, the changes in diffusivity vary significantly among the different simulations. For example, simulations with subsurface canopies show a large increase in diffusivity in the upper portion of the farm (both within and between the rows; Supplementary Figures 2, 4, 6), likely produced by the strong shear layer that develops at the top of the farm in the absence of a floating dense canopy.

In general, there is agreement between kelp simulated by the MAG model and field observations of M. pyrifera grown on longline and from natural kelp forests in close proximity to the simulated farm. The average growth rate of MAG is 4% d^-1 (Figure 3D). Also, the biomass per fronds, the distribution of biomass among fronds, the total number of fronds per plant, and the total biomass per plant agree with data from plants sampled by the Santa Barbara Coastal Long Term Ecological Research program (Figures 3B,C).

To simulate the physical effects of farm drag on the nearshore environment, we add an approximation of a vertically variable kelp-drag to ROMS. The drag force due to kelp is parameterized with Equation (5). In ROMS, F_D acts as a sink on horizontal momentum.

For illustrative purposes, we choose a single frond surface area per unit volume, a(z), profile that represents a maximal kelp canopy, where frond surface area density is largest near the surface (Supplementary Table 3). For this demonstration, we also utilize a higher-resolution (Δx = 36 m) ROMS hindcast simulation of the Santa Barbara Channel without biogeochemistry, described in Dauhajre et al. (2019), that better resolves both smaller-scale, turbulent nearshore currents and farm-flow interactions as ambient currents flow through and around the farm. In this simulation, we define a 432 x 432 m farm to match the resolution of the ROMS hindcast with homogenous a(z) within the farm (Figure 2A). The hindcast spans a period of 10 days in November 2006. Continual release of a passive tracer within the farm boundaries demonstrates the effects of farm drag on local shelf flows and calculation of a “footprint” of farm effluent over the 10-day period.

Figure 10 demonstrates the effects of the (parameterized) farm drag with visualization of a passive tracer with (unitless) concentration C that is continuously released within the farm. In this manner, the tracer acts as a stand-in for any farm “effluent.” Such effluent will be advected by the shelf currents (altered by the farm itself) and may reach the shoreline or influence local ecosystem functioning. Physically, the farm drag dampens currents and causes a diversion of incoming flow beneath the farm (not shown). Additionally, the interaction of farm drag and ambient shelf currents often generates wake vortices (Figure 10A) that form downstream of the mean flow and live for hours to days. These vortices are analogous to headland wakes, exhibit large lateral shear in velocity, and influence the redistribution of farm effluent on the shelf. The snapshot in Figure 10A illustrates how these “farm wakes” advect the vertically averaged farm effluent (〈C〉_z) which can reach ~ 10–15 km downstream of the farm.

A time-averaged view of the farm effluent distribution shows a substantial farm “footprint” due to a primarily westward, along-shore flow during the 10-day simulation period. We plot this footprint as the temporal root-mean square of the vertically integrated tracer concentration in the upper 20 m (σ_〈C〉z, where 〈.〉 denotes a 10-day temporal root-mean square), shown in Figure 10B. Generally, the farm effluent extends upcoast (westward) of the farm, due to the predominance of a large-scale, westward along-shore flow. However, as the snapshot (Figure 10A) illustrates, more complex, smaller-scale flow patterns ultimately contribute to the mean distribution of farm effluent. Notably, 10% of the total released effluent (solid gray line, Figure 10B) extends ~ 20 km in the along-shore direction. This relatively large footprint suggests non-trivial, significant affects of the farm on the continental shelf environment due to a combination of local flow alteration by drag and associated transport of farm by-products on the shelf. While we limit the analysis here to a fully grown stage farm with a large canopy on the surface, we observe similar footprints for other a(z) profiles with reduced surface canopy. Similarly, we leave more detailed investigation of local ecosystem effects, with a coupled ROMS-BEC-MAG model, for future study.