diff --git a/src/Fields/broadcast.jl b/src/Fields/broadcast.jl
index 14cb646eb4..a63c174e2a 100644
--- a/src/Fields/broadcast.jl
+++ b/src/Fields/broadcast.jl
@@ -294,6 +294,13 @@ Base.Broadcast.broadcasted(
 ) where {T <: Geometry.AxisTensor} =
     Base.Broadcast.broadcasted(fs, (x...) -> T(x...), args...)
 
+Base.Broadcast.broadcasted(
+    ::typeof(Base.literal_pow),
+    ::typeof(^),
+    ::Field,
+    ::Val{n},
+) where {n} = Base.Broadcast.broadcasted(x -> Base.literal_pow(^, x, Val(n)), f)
+
 # Specialize handling of +, *, muladd, so that we can support broadcasting over NamedTuple element types
 # Required for ODE solvers
 
diff --git a/test/Fields/field.jl b/test/Fields/field.jl
index ac1f730295..673054f64a 100644
--- a/test/Fields/field.jl
+++ b/test/Fields/field.jl
@@ -94,6 +94,21 @@ end
     @test axes(point_field) isa Spaces.PointSpace
 end
 
+# https://github.com/CliMA/ClimaCore.jl/issues/1126
+function pow_n(f)
+    @. f.x = f.x^2
+    return nothing
+end
+@testset "Broadcasting with ^n" begin
+    FT = Float32
+    for space in TU.all_spaces(FT)
+        f = fill((; x = FT(1)), space)
+        pow_n(f) # Compile first
+        p_allocated = @allocated pow_n(f)
+        @test p_allocated == 0
+    end
+end
+
 # Requires `--check-bounds=yes`
 @testset "Constructing & broadcasting over empty fields" begin
     FT = Float32