diff --git a/Project.toml b/Project.toml
index ab8caccf4..1ed7a18d2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "KernelFunctions"
 uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
-version = "0.10.20"
+version = "0.10.21"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
diff --git a/docs/src/kernels.md b/docs/src/kernels.md
index 7734831f2..56f3cfb9b 100644
--- a/docs/src/kernels.md
+++ b/docs/src/kernels.md
@@ -27,6 +27,7 @@ CosineKernel
 
 ```@docs
 ExponentialKernel
+GibbsKernel
 LaplacianKernel
 SqExponentialKernel
 SEKernel
diff --git a/src/KernelFunctions.jl b/src/KernelFunctions.jl
index c0290b390..28271cb47 100644
--- a/src/KernelFunctions.jl
+++ b/src/KernelFunctions.jl
@@ -17,6 +17,7 @@ export PiecewisePolynomialKernel
 export PeriodicKernel, NeuralNetworkKernel
 export KernelSum, KernelProduct, KernelTensorProduct
 export TransformedKernel, ScaledKernel, NormalizedKernel
+export GibbsKernel
 
 export Transform,
     SelectTransform,
@@ -53,7 +54,7 @@ using Functors
 using LinearAlgebra
 using Requires
 using SpecialFunctions: loggamma, besselk, polygamma
-using IrrationalConstants: logtwo, twoπ
+using IrrationalConstants: logtwo, twoπ, invsqrt2
 using LogExpFunctions: softplus
 using StatsBase
 using TensorCore
@@ -94,6 +95,7 @@ include("basekernels/rational.jl")
 include("basekernels/sm.jl")
 include("basekernels/wiener.jl")
 
+include("kernels/gibbskernel.jl")
 include("kernels/scaledkernel.jl")
 include("kernels/normalizedkernel.jl")
 include("matrix/kernelmatrix.jl")
diff --git a/src/kernels/gibbskernel.jl b/src/kernels/gibbskernel.jl
new file mode 100644
index 000000000..4ce416f7f
--- /dev/null
+++ b/src/kernels/gibbskernel.jl
@@ -0,0 +1,42 @@
+"""
+    GibbsKernel(; lengthscale)
+
+# Definition
+
+The Gibbs kernel is non-stationary generalisation of the squared exponential
+kernel. The lengthscale parameter ``l`` becomes a function of
+position ``l(x)``.
+
+For a constant function``l(x) = c``, one recovers the standard squared exponential kernel
+with lengthscale `c`.
+
+```math
+k(x, y; l) = \\sqrt{ \\left(\\frac{2 l(x) l(y)}{l(x)^2 + l(y)^2} \\right) }
+\\quad \\rm{exp} \\left( - \\frac{(x - y)^2}{l(x)^2 + l(y)^2} \\right)
+```
+
+# References
+
+Mark N. Gibbs. "Bayesian Gaussian Processes for Regression and Classication." PhD thesis, 1997
+
+Christopher J. Paciorek and Mark J. Schervish. "Nonstationary Covariance Functions
+for Gaussian Process Regression". NeurIPS, 2003
+
+Sami Remes, Markus Heinonen, Samuel Kaski. "Non-Stationary Spectral Kernels". arXiV:1705.08736, 2017
+
+Sami Remes, Markus Heinonen, Samuel Kaski. "Neural Non-Stationary Spectral Kernel". arXiv:1811.10978, 2018
+"""
+struct GibbsKernel{T} <: Kernel
+    lengthscale::T
+end
+
+GibbsKernel(; lengthscale) = GibbsKernel(lengthscale)
+
+function (k::GibbsKernel)(x, y)
+    lengthscale = k.lengthscale
+    lx = lengthscale(x)
+    ly = lengthscale(y)
+    l = invsqrt2 * hypot(lx, ly)
+    kernel = (sqrt(lx * ly) / l) * with_lengthscale(SqExponentialKernel(), l)
+    return kernel(x, y)
+end
diff --git a/test/kernels/gibbskernel.jl b/test/kernels/gibbskernel.jl
new file mode 100644
index 000000000..3c1722dcd
--- /dev/null
+++ b/test/kernels/gibbskernel.jl
@@ -0,0 +1,13 @@
+@testset "gibbskernel" begin
+    x = randn()
+    y = randn()
+
+    # this is the gibbs lengthscale function.
+    ell(x) = exp(sum(sin, x))
+    # create a gibbs kernel with our specific lengthscale function
+    k_gibbs = GibbsKernel(ell)
+
+    @test k_gibbs(x, y) ≈
+          sqrt((2 * ell(x) * ell(y)) / (ell(x)^2 + ell(y)^2)) *
+          exp(-(x - y)^2 / (ell(x)^2 + ell(y)^2))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 2ace8981f..63c8e8f89 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -126,6 +126,7 @@ include("test_utils.jl")
             include("kernels/transformedkernel.jl")
             include("kernels/normalizedkernel.jl")
             include("kernels/neuralkernelnetwork.jl")
+            include("kernels/gibbskernel.jl")
         end
         @info "Ran tests on Kernel"
     end