diff --git a/examples/global_climate_simulation.jl b/examples/global_climate_simulation.jl
index 0f461e63..b3bdc66b 100644
--- a/examples/global_climate_simulation.jl
+++ b/examples/global_climate_simulation.jl
@@ -49,7 +49,7 @@ tracer_advection     = WENO(order=5)
 free_surface         = SplitExplicitFreeSurface(grid; substeps=40)
 catke_closure        = NumericalEarth.Oceans.default_ocean_closure()
 eddy_closure         = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3)
-viscous_closure      = Oceananigans.TurbulenceClosures.HorizontalScalarBiharmonicDiffusivity(ν=1e12)
+viscous_closure      = area_scaled_biharmonic_viscosity()
 closures             = (catke_closure, eddy_closure, viscous_closure)
 nothing #hide
 
diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl
index 1e50a243..dc2ac5e8 100644
--- a/examples/one_degree_simulation.jl
+++ b/examples/one_degree_simulation.jl
@@ -48,12 +48,13 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height);
 # ### Closures
 #
 # We include a Gent-McWilliams isopycnal diffusivity as a parameterization for the mesoscale
-# eddy fluxes. For vertical mixing at the upper-ocean boundary layer we include the CATKE
-# parameterization.
+# eddy fluxes, an area-scaled biharmonic viscosity for grid-scale dissipation, and the CATKE
+# parameterization for vertical mixing at the upper-ocean boundary layer.
 
 using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity, AdvectiveFormulation
 
 eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3, skew_flux_formulation=AdvectiveFormulation())
+horizontal_viscosity = area_scaled_biharmonic_viscosity()
 vertical_mixing = NumericalEarth.Oceans.default_ocean_closure()
 
 # ### Ocean simulation
@@ -65,7 +66,7 @@ momentum_advection = WENOVectorInvariant(order=5)
 tracer_advection   = WENO(order=5)
 
 ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface,
-                         closure=(eddy_closure, vertical_mixing))
+                         closure=(eddy_closure, horizontal_viscosity, vertical_mixing))
 
 @info "We've built an ocean simulation with model:"
 @show ocean.model
diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl
index 2ef5330c..1ec03db5 100644
--- a/src/NumericalEarth.jl
+++ b/src/NumericalEarth.jl
@@ -89,6 +89,7 @@ export
     DatasetRestoring,
     atmosphere_simulation,
     ocean_simulation,
+    area_scaled_biharmonic_viscosity,
     sea_ice_simulation,
     default_sea_ice,
     sea_ice_dynamics,
diff --git a/src/Oceans/Oceans.jl b/src/Oceans/Oceans.jl
index 7e31455c..e8470d95 100644
--- a/src/Oceans/Oceans.jl
+++ b/src/Oceans/Oceans.jl
@@ -1,6 +1,6 @@
 module Oceans
 
-export ocean_simulation, SlabOcean
+export ocean_simulation, SlabOcean, area_scaled_biharmonic_viscosity
 
 using Adapt: Adapt, adapt
 using KernelAbstractions: @kernel, @index
@@ -16,12 +16,12 @@ using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, ImmersedBoundaryCon
 using Oceananigans.Models.HydrostaticFreeSurfaceModels: HydrostaticFreeSurfaceModel
 using Oceananigans.Models.HydrostaticFreeSurfaceModels.SplitExplicitFreeSurfaces: SplitExplicitFreeSurface
 using Oceananigans.OrthogonalSphericalShellGrids: OrthogonalSphericalShellGrids, TripolarGrid
-using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ
+using Oceananigans.Operators: ℑxyᶠᶜᵃ, ℑxyᶜᶠᵃ, Az
 using Oceananigans.Simulations: Simulation
 using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity,
                                                                      CATKEMixingLength,
                                                                      CATKEEquation
-using Oceananigans.Units: minutes, hours
+using Oceananigans.Units: minutes, hours, days
 using Oceananigans.Utils: with_tracers, launch!
 using SeawaterPolynomials: SeawaterPolynomials
 using SeawaterPolynomials.TEOS10: TEOS10EquationOfState
diff --git a/src/Oceans/ocean_simulation.jl b/src/Oceans/ocean_simulation.jl
index dbb2453f..31b64d84 100644
--- a/src/Oceans/ocean_simulation.jl
+++ b/src/Oceans/ocean_simulation.jl
@@ -3,7 +3,8 @@ using Oceananigans.BoundaryConditions: DefaultBoundaryCondition
 using Oceananigans.DistributedComputations: DistributedGrid, all_reduce
 using Oceananigans.Grids: inactive_node
 using Oceananigans.OrthogonalSphericalShellGrids
-using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization
+using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization,
+                                       HorizontalScalarBiharmonicDiffusivity
 using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity,
                                                                      CATKEMixingLength,
                                                                      CATKEEquation
@@ -84,6 +85,24 @@ function default_ocean_closure(FT=Oceananigans.defaults.FloatType)
     return CATKEVerticalDiffusivity(VerticallyImplicitTimeDiscretization(), FT; mixing_length, turbulent_kinetic_energy_equation)
 end
 
+@inline _area_scaled_biharmonic_viscosity(i, j, k, grid, ℓx, ℓy, ℓz, clock, fields, λ) =
+    Az(i, j, k, grid, ℓx, ℓy, ℓz)^2 / λ
+
+"""
+    area_scaled_biharmonic_viscosity(FT=Oceananigans.defaults.FloatType; timescale=15days)
+
+Return a `HorizontalScalarBiharmonicDiffusivity` with a grid-area-scaled viscosity
+coefficient `ν = Az² / timescale`, where `Az` is the horizontal area of each grid cell.
+This produces a viscosity that adapts to local grid cell size, providing appropriate
+damping at all latitudes.
+"""
+function area_scaled_biharmonic_viscosity(FT=Oceananigans.defaults.FloatType; timescale=15days)
+    return HorizontalScalarBiharmonicDiffusivity(FT;
+        ν = _area_scaled_biharmonic_viscosity,
+        discrete_form = true,
+        parameters = timescale)
+end
+
 # Two-band shortwave penetration in the Paulson & Simpson (1977) form,
 # Defaults are Jerlov Type I (clearest open-ocean water)
 function default_radiative_forcing(grid)