diff --git a/src/quantities.jl b/src/quantities.jl
index 1024158c..64139c63 100644
--- a/src/quantities.jl
+++ b/src/quantities.jl
@@ -265,17 +265,17 @@ isapprox(x::Number, y::AbstractQuantity; kwargs...) = isapprox(y, x; kwargs...)
 function isapprox(
     x::AbstractArray{<:AbstractQuantity{T1,D,U1}},
     y::AbstractArray{<:AbstractQuantity{T2,D,U2}};
-    rtol::Real=Base.rtoldefault(T1,T2,0),
     atol=zero(Quantity{real(T1),D,U1}),
+    rtol::Real=Base.rtoldefault(T1,T2,atol>zero(atol)),
+    nans::Bool=false,
     norm::Function=norm,
 ) where {T1,D,U1,T2,U2}
-
     d = norm(x - y)
     if isfinite(d)
-        return d <= atol + rtol*max(norm(x), norm(y))
+        return iszero(rtol) ? d <= atol : d <= max(atol, rtol*max(norm(x), norm(y)))
     else
         # Fall back to a component-wise approximate comparison
-        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol), zip(x, y))
+        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol, nans=nans), zip(x, y))
     end
 end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index b7a2e39c..e97e6529 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1542,6 +1542,10 @@ end
             @test isapprox([1cm, 200cm], [0.01m, 2.0m])
             @test !isapprox([1.0], [1.0m])
             @test !isapprox([1.0m], [1.0])
+            @test isapprox([1.0m, NaN*m], [nextfloat(1.0)*m, NaN*m], nans=true)
+            @test !isapprox([1.0m, NaN*m], [nextfloat(1.0)*m, NaN*m], nans=false)
+            @test !isapprox([1.0m, 2.0m], [1.1m, 2.2m], rtol=0.05, atol=0.2m)
+            @test !isapprox([1.0m], [nextfloat(1.0)*m], atol=eps(0.1)*m)
         end
         @testset ">> Unit stripping" begin
             @test @inferred(ustrip([1u"m", 2u"m"])) == [1,2]