NumericalEarth · simone-silvestri · Apr 16, 2026 · Apr 16, 2026 · Apr 16, 2026 · Apr 16, 2026
diff --git a/Project.toml b/Project.toml
@@ -13,6 +13,8 @@ DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
 Glob = "c27321d9-0574-5035-807b-f59d2c89b15c"
@@ -65,6 +67,8 @@ DataDeps = "0.7"
 Dates = "<0.0.1, 1"
 Distances = "0.10"
 DocStringExtensions = "0.9"
+Documenter = "1.17.0"
+DocumenterCitations = "1.4.1"
 Downloads = "<0.0.1, 1"
 GPUArraysCore = "0.2.0"
 Glob = "1"

diff --git a/docs/src/earth_system_model.md b/docs/src/earth_system_model.md
@@ -120,7 +120,7 @@ To change the defaults, construct `ComponentInterfaces` yourself and pass it in.
 For example, to use constant transfer coefficients instead of similarity theory:
 
 ```jldoctest esm
-atmosphere_ocean_fluxes = CoefficientBasedFluxes(drag_coefficient=2e-3)
+atmosphere_ocean_fluxes = CoefficientBasedFluxes(transfer_coefficients = (2e-3, 2e-3, 2e-3))
 interfaces = NumericalEarth.EarthSystemModels.ComponentInterfaces(nothing, ocean;
                                                                   atmosphere_ocean_fluxes)
 model = OceanOnlyModel(ocean; interfaces)

diff --git a/docs/src/interface_fluxes.md b/docs/src/interface_fluxes.md
@@ -137,9 +137,31 @@ A comprehensive example is given below, but we note briefly here that
 ```@example interface_fluxes
 using NumericalEarth
 
-coefficient_fluxes = CoefficientBasedFluxes(drag_coefficient=2e-3,
-                                            heat_transfer_coefficient=2e-3,
-                                            vapor_flux_coefficient=1e-3)
+coefficient_fluxes = CoefficientBasedFluxes(transfer_coefficients = (2e-3, 2e-3, 1e-3))
+```
+
+Alternatively, the drag coefficient can be specified as a wind-speed-dependent polynomial
+following Large & Yeager (2004). In this case `CoefficientBasedFluxes` evaluates the polynomial
+at each iteration rather than using a constant:
+
+```@example interface_fluxes
+using NumericalEarth.EarthSystemModels.InterfaceComputations: PolynomialNeutralDragCoefficient
+
+poly_drag = PolynomialNeutralDragCoefficient()
+poly_fluxes = CoefficientBasedFluxes(transfer_coefficients = (poly_drag, 1e-3, 1e-3))
+```
+
+For the full Large & Yeager (2004) bulk algorithm with stability corrections,
+use `LargeYeagerTransferCoefficients`. This computes all three transfer coefficients
+(drag, Stanton, Dalton) from the neutral drag polynomial with Monin-Obukhov
+stability corrections (L&Y eqs. 6c-6d, 10a-10c):
+
+```@example interface_fluxes
+using NumericalEarth.EarthSystemModels.InterfaceComputations: FixedIterations, LargeYeagerTransferCoefficients
+
+ly = LargeYeagerTransferCoefficients()
+ly_fluxes = CoefficientBasedFluxes(transfer_coefficients = ly,
+                                   solver_stop_criteria = FixedIterations(5))
 ```
 
 ### Similarity theory for neutral boundary layers
@@ -335,13 +357,18 @@ neutral_similarity_fluxes = SimilarityTheoryFluxes(stability_functions=nothing;
 interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=neutral_similarity_fluxes)
 increased_roughness_model = OceanOnlyModel(ocean; atmosphere, interfaces)
 
-coefficient_fluxes = CoefficientBasedFluxes(drag_coefficient=2e-3)
+coefficient_fluxes = CoefficientBasedFluxes(transfer_coefficients = (2e-3, 2e-3, 2e-3))
 interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=coefficient_fluxes)
 coefficient_model = OceanOnlyModel(ocean; atmosphere, interfaces)
+
+ly_fluxes = CoefficientBasedFluxes(transfer_coefficients = LargeYeagerTransferCoefficients(),
+                                   solver_stop_criteria = FixedIterations(5))
+interfaces = ComponentInterfaces(atmosphere, ocean; atmosphere_ocean_fluxes=ly_fluxes)
+ly_model = OceanOnlyModel(ocean; atmosphere, interfaces)
 ```
 
 Note that `EarthSystemModel` computes fluxes upon instantiation, so after constructing
-the two models we are ready to analyze the results.
+the models we are ready to analyze the results.
 We first verify that the similarity model friction velocity has been computed successfully,
 
 ```@example interface_fluxes
@@ -370,7 +397,16 @@ Cᴰ_coeff = @. (u★_coeff / uᵃᵗ)^2
 extrema(Cᴰ_coeff)
 ```
 
-We'll compare the computed fluxes and drag coefficients from our two models with
+We also extract the drag coefficient from the Large & Yeager transfer coefficient model:
+
+```@example interface_fluxes
+u★_ly = ly_model.interfaces.atmosphere_ocean_interface.fluxes.friction_velocity
+u★_ly = interior(u★_ly, :, 1, 1)
+Cᴰ_ly = @. (u★_ly / uᵃᵗ)^2
+extrema(Cᴰ_ly)
+```
+
+We'll compare the computed fluxes and drag coefficients from our models with
 a polynomial expression due to [large2009global](@citet), and
 an expression reported by [edson2013exchange](@citet) that was developed at ECMWF,
 
@@ -406,16 +442,18 @@ Cᴰ_rough = @. (u★_rough / uᵃᵗ)^2
 
 fig = Figure(size=(800, 400))
 axu = Axis(fig[1:2, 1], xlabel="uᵃᵗ (m s⁻¹) at 10 m", ylabel="u★ (m s⁻¹)")
-lines!(axu, uᵃᵗ, u★, label="NumericalEarth default")
-lines!(axu, uᵃᵗ, u★_rough, label="Increased roughness model")
-lines!(axu, uᵃᵗ, u★_LY, label="Large and Yeager (2009) polynomial fit")
-lines!(axu, uᵃᵗ, u★_EC, label="ECMWF polynomial fit (Edson et al. 2013)")
+lines!(axu, uᵃᵗ, u★, label="SimilarityTheoryFluxes (default)")
+lines!(axu, uᵃᵗ, u★_rough, label="SimilarityTheoryFluxes (increased roughness)")
+lines!(axu, uᵃᵗ, u★_ly, label="LargeYeagerTransferCoefficients (L&Y 2004)")
+lines!(axu, uᵃᵗ, u★_LY, label="Large and Yeager (2009) polynomial fit", linestyle=:dash)
+lines!(axu, uᵃᵗ, u★_EC, label="ECMWF polynomial fit (Edson et al. 2013)", linestyle=:dash)
 
 axd = Axis(fig[1:2, 2], xlabel="uᵃᵗ (m s⁻¹) at 10 m", ylabel="1000 × Cᴰ")
-lines!(axd, uᵃᵗ, 1000 .* Cᴰ_default, label="NumericalEarth default")
-lines!(axd, uᵃᵗ, 1000 .* Cᴰ_rough, label="Increased roughness model")
-lines!(axd, uᵃᵗ, 1000 .* Cᴰ_LY, label="Large and Yeager (2009) polynomial fit")
-lines!(axd, uᵃᵗ, 1000 .* Cᴰ_EC, label="ECMWF polynomial fit (Edson et al. 2013)")
+lines!(axd, uᵃᵗ, 1000 .* Cᴰ_default, label="SimilarityTheoryFluxes (default)")
+lines!(axd, uᵃᵗ, 1000 .* Cᴰ_rough, label="SimilarityTheoryFluxes (increased roughness)")
+lines!(axd, uᵃᵗ, 1000 .* Cᴰ_ly, label="LargeYeagerTransferCoefficients (L&Y 2004)")
+lines!(axd, uᵃᵗ, 1000 .* Cᴰ_LY, label="Large and Yeager (2009) polynomial fit", linestyle=:dash)
+lines!(axd, uᵃᵗ, 1000 .* Cᴰ_EC, label="ECMWF polynomial fit (Edson et al. 2013)", linestyle=:dash)
 
 Legend(fig[3, 1:2], axd, nbanks = 2)
 

diff --git a/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl b/src/EarthSystemModels/InterfaceComputations/InterfaceComputations.jl
@@ -13,10 +13,15 @@ export
     MomentumRoughnessLength,
     ScalarRoughnessLength,
     CoefficientBasedFluxes,
+    SimilarityScales,
+    PolynomialNeutralDragCoefficient,
+    LargeYeagerTransferCoefficients,
+    LinearStableStabilityFunction,
     SkinTemperature,
     BulkTemperature,
     atmosphere_ocean_stability_functions,
     atmosphere_sea_ice_stability_functions,
+    large_yeager_stability_functions,
     compute_atmosphere_ocean_fluxes!,
     compute_atmosphere_sea_ice_fluxes!,
     compute_sea_ice_ocean_fluxes!,