diff --git a/src/covariance.jl b/src/covariance.jl
index dcb84fe..a4bb283 100644
--- a/src/covariance.jl
+++ b/src/covariance.jl
@@ -77,6 +77,23 @@ Evaluate covariance `cov` between all elements of `domain₁` and `domain₂`.
 """
 pairwise(cov::Covariance, args...) = sill(cov.γ) .- pairwise(cov.γ, args...)
 
+"""
+    pairwise!(C, cov, domain)
+
+Evaluates covariance `cov` between all elements in the `domain` in-place, filling the matrix `C`.
+
+    pairwise!(C, cov, domain₁, domain₂)
+
+Evaluates covariance `cov` between all elements of `domain₁` and `domain₂` in-place, filling the matrix `C`.
+"""
+function pairwise!(C, cov::Covariance, args...)
+  pairwise!(C, cov.γ, args...)
+  for i in eachindex(Γ)
+    C[i] = sill(cov.γ) - C[i]
+  end
+  C
+end
+
 # -----------
 # IO METHODS
 # -----------
diff --git a/src/variogram.jl b/src/variogram.jl
index 5818194..de09bec 100644
--- a/src/variogram.jl
+++ b/src/variogram.jl
@@ -152,6 +152,11 @@ function pairwise(γ::Variogram, domain)
   pairwise!(Γ, γ, domain)
 end
 
+"""
+    pairwise!(Γ, γ, domain)
+
+Evaluates covariance `γ` between all elements in the `domain` in-place, filling the matrix `Γ`.
+"""
 function pairwise!(Γ, γ::Variogram, domain)
   n = length(domain)
   @inbounds for j in 1:n
@@ -185,6 +190,11 @@ function pairwise(γ::Variogram, domain₁, domain₂)
   pairwise!(Γ, γ, domain₁, domain₂)
 end
 
+"""
+    pairwise!(Γ, γ, domain₁, domain₂)
+
+Evaluates covariance `γ` between all elements of `domain₁` and `domain₂` in-place, filling the matrix `Γ`.
+"""
 function pairwise!(Γ, γ::Variogram, domain₁, domain₂)
   m = length(domain₁)
   n = length(domain₂)