diff --git a/NEWS.md b/NEWS.md
index 1a50ef7..05a5824 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -253,6 +253,27 @@ julia> units(se_var)
 "K^2"
 ```
 
+### Plotting bias
+
+To plot the bias and display information such as the root mean squared error (RMSE) and the
+global bias in the plot, you use the function `plot_bias!(fig, sim, obs)`. In the example
+below, we plot the bias between our simulation and some observations stored in
+"ta\_1d\_average.nc".
+
+```julia
+import ClimaAnalysis
+import ClimaAnalysis.Visualize: plot_bias!
+import GeoMakie
+import CairoMakie
+
+obs_var = ClimaAnalysis.OutputVar("ta_1d_average.nc")
+sim_var = ClimaAnalysis.get(ClimaAnalysis.simdir("simulation_output"), "ta")
+
+fig = CairoMakie.Figure()
+plot_bias!(fig, sim_var, obs_var)
+CairoMakie.save("myfigure.pdf", fig)
+```
+
 ## Bug fixes
 
 - Increased the default value for `warp_string` to 72.
diff --git a/docs/src/api.md b/docs/src/api.md
index b39d372..684a2f3 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -106,4 +106,5 @@ Visualize.oceanmask
 Visualize.landmask
 Visualize.contour2D_on_globe!
 Visualize.heatmap2D_on_globe!
+Visualize.plot_bias!
 ```
diff --git a/docs/src/assets/bias_plot.png b/docs/src/assets/bias_plot.png
new file mode 100644
index 0000000..d09ec70
Binary files /dev/null and b/docs/src/assets/bias_plot.png differ
diff --git a/docs/src/visualize.md b/docs/src/visualize.md
index 9e62fa8..4593485 100644
--- a/docs/src/visualize.md
+++ b/docs/src/visualize.md
@@ -45,3 +45,29 @@ more easily.
 The output might look something like:
 
 ![oceanmask](./assets/oceanmask.png)
+
+### Plotting bias
+
+After [computing the bias](@ref bias) between observational and simulation data, you may
+want to plot the bias and display information such as the root mean squared error (RMSE) and
+the global bias in the plot. To do this, you use the function [`plot_bias!(fig, sim,
+obs)`](@ref Visualize.plot_bias!). In the example below, we plot the bias between our
+simulation and some observations stored in "ta\_1d\_average.nc".
+
+```julia
+import ClimaAnalysis
+import ClimaAnalysis.Visualize: plot_bias!
+import GeoMakie
+import CairoMakie
+
+obs_var = ClimaAnalysis.OutputVar("ta_1d_average.nc")
+sim_var = ClimaAnalysis.get(ClimaAnalysis.simdir("simulation_output"), "ta")
+
+fig = CairoMakie.Figure()
+plot_bias!(fig, sim_var, obs_var)
+CairoMakie.save("myfigure.pdf", fig)
+```
+
+The output produces something like:
+
+![biasplot](./assets/bias_plot.png)
\ No newline at end of file
diff --git a/ext/ClimaAnalysisGeoMakieExt.jl b/ext/ClimaAnalysisGeoMakieExt.jl
index 4b74514..7c108dc 100644
--- a/ext/ClimaAnalysisGeoMakieExt.jl
+++ b/ext/ClimaAnalysisGeoMakieExt.jl
@@ -254,4 +254,125 @@ function Visualize.contour2D_on_globe!(
     )
 end
 
+"""
+    plot_bias!(fig::Makie.Figure,
+               sim::ClimaAnalysis.OutputVar,
+               obs::ClimaAnalysis.OutputVar;
+               cmap_extrema = extrema(ClimaAnalysis.bias(sim, obs).data),
+               p_loc = (1, 1),
+               plot_coastline = true,
+               plot_colorbar = true,
+               mask = nothing,
+               more_kwargs)
+    plot_bias!(grid_layout::Makie.GridLayout,
+               sim::ClimaAnalysis.OutputVar,
+               obs::ClimaAnalysis.OutputVar;
+               cmap_extrema = extrema(ClimaAnalysis.bias(sim, obs).data),
+               p_loc = (1, 1),
+               plot_coastline = true,
+               plot_colorbar = true,
+               mask = nothing,
+               more_kwargs)
+
+
+Plot the bias (`sim.data - var.data`) on a projected geoid. The gloal bias and root mean
+squared error (RMSE) are computed and can be found in the title of the plot. This function
+plots the returned `OutputVar` of `ClimaAnalysis.bias(sim, obs)`. See also
+[`ClimaAnalysis.bias`](@ref).
+
+The plot comes with labels, units, and a colorbar. This function uses a constrained colormap
+based on the values of `cmap_extrema`.
+
+The dimensions have to be longitude and latitude.
+
+`mask` has to be an object that can be plotted by `Makie.poly`. `ClimaAnalysis` comes with
+predefined masks, check out [`Visualize.oceanmask`](@ref) and [`Visualize.landmask`](@ref).
+
+!!! note
+    Masking does not affect the colorbar. If you have values defined beneath the map, they
+    can still affect the colorbar.
+
+Additional arguments to the plotting and axis functions
+=======================================================
+
+`more_kwargs` can be a dictionary that maps symbols to additional options for:
+- the axis (`:axis`)
+- the plotting function (`:plot`)
+- the colorbar (`:cb`)
+- the coastline (`:coast`)
+- the mask (`:mask`)
+
+The coastline is plotted from `GeoMakie.coastline` using the `lines!` plotting function.
+
+The values are splatted in the relevant functions. Populate them with a
+Dictionary of `Symbol`s => values to pass additional options.
+"""
+function Visualize.plot_bias!(
+    place::MakiePlace,
+    sim::ClimaAnalysis.OutputVar,
+    obs::ClimaAnalysis.OutputVar;
+    cmap_extrema = extrema(ClimaAnalysis.bias(sim, obs).data),
+    p_loc = (1, 1),
+    plot_coastline = true,
+    plot_colorbar = true,
+    mask = nothing,
+    more_kwargs = Dict(
+        :plot => Dict(),
+        :cb => Dict(),
+        :axis => Dict(),
+        :coast => Dict(:color => :black),
+        :mask => Dict(),
+    ),
+)
+    bias_var = ClimaAnalysis.bias(sim, obs)
+    global_bias = round(bias_var.attributes["global_bias"], sigdigits = 3)
+    rmse = round(ClimaAnalysis.global_rmse(sim, obs), sigdigits = 3)
+    units = ClimaAnalysis.units(bias_var)
+
+    bias_var.attributes["long_name"] *= " (RMSE: $rmse $units, Global bias: $global_bias $units)"
+    min_level, max_level = cmap_extrema
+
+    # Make sure that 0 is at the center
+    cmap = Visualize._constrained_cmap(
+        Makie.cgrad(:vik).colors,
+        min_level,
+        max_level;
+        categorical = true,
+    )
+    nlevels = 11
+    # Offset so that it covers 0
+    levels = collect(range(min_level, max_level, length = nlevels))
+    offset = levels[argmin(abs.(levels))]
+    levels = levels .- offset
+    ticklabels = map(x -> string(round(x; digits = 0)), levels)
+    ticks = (levels, ticklabels)
+
+    default_kwargs = Dict(
+        :plot => Dict(
+            :colormap => cmap,
+            :levels => levels,
+            :extendhigh => :auto,
+            :extendlow => :auto,
+        ),
+        :cb => Dict(:ticks => ticks),
+    )
+
+    # Function for recursively merging two dictionaries if the values of the dictionaries
+    # are dictionaries and the values of those are also dictionaries and so on
+    # See: https://discourse.julialang.org/t/multi-layer-dict-merge/27261/6
+    recursive_merge(x::AbstractDict...) = merge(recursive_merge, x...)
+    recursive_merge(x...) = x[end]
+    default_and_more_kwargs = recursive_merge(default_kwargs, more_kwargs)
+
+    return Visualize.contour2D_on_globe!(
+        place,
+        bias_var;
+        p_loc = p_loc,
+        plot_coastline = plot_coastline,
+        plot_colorbar = plot_colorbar,
+        mask = mask,
+        more_kwargs = default_and_more_kwargs,
+    )
+end
+
 end
diff --git a/src/Visualize.jl b/src/Visualize.jl
index e2e59ae..ac770bd 100644
--- a/src/Visualize.jl
+++ b/src/Visualize.jl
@@ -28,4 +28,6 @@ function contour2D_on_globe! end
 
 function heatmap2D_on_globe! end
 
+function plot_bias! end
+
 end
diff --git a/test/test_GeoMakieExt.jl b/test/test_GeoMakieExt.jl
index 5f5c49c..ce5f3a8 100644
--- a/test/test_GeoMakieExt.jl
+++ b/test/test_GeoMakieExt.jl
@@ -100,4 +100,40 @@ using OrderedCollections
 
     output_name = joinpath(tmp_dir, "test_contours2D_globe_with_oceanmask.png")
     Makie.save(output_name, fig5)
+
+    # Test plot_bias
+    fig6 = Makie.Figure()
+
+    lon = collect(range(-179.5, 179.5, 360))
+    lat = collect(range(-89.5, 89.5, 180))
+    data = collect(reshape(-32400:32399, (360, 180))) ./ (32399.0 / 5.0)
+    dims = OrderedDict(["lon" => lon, "lat" => lat])
+    attribs =
+        Dict("long_name" => "idk", "short_name" => "short", "units" => "kg")
+    dim_attribs = OrderedDict([
+        "lon" => Dict("units" => "deg"),
+        "lat" => Dict("units" => "deg"),
+    ])
+    var = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data)
+
+    data_zero = zeros(length(lon), length(lat))
+    var_zero = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, data_zero)
+
+    ClimaAnalysis.Visualize.plot_bias!(fig6, var, var_zero)
+    output_name = joinpath(tmp_dir, "plot_bias.png")
+    Makie.save(output_name, fig6)
+
+    # Test plot bias with keyword arguments
+    fig7 = Makie.Figure()
+    ClimaAnalysis.Visualize.plot_bias!(
+        fig7,
+        var,
+        var_zero,
+        more_kwargs = Dict(
+            :axis => Dict(:title => "no title"),
+            :plot => Dict(:extendhigh => nothing),
+        ),
+    )
+    output_name = joinpath(tmp_dir, "plot_bias_kwargs.png")
+    Makie.save(output_name, fig7)
 end