diff --git a/act/plotting/contourdisplay.py b/act/plotting/contourdisplay.py
index cbfe7ef69d..75e8bebdd3 100644
--- a/act/plotting/contourdisplay.py
+++ b/act/plotting/contourdisplay.py
@@ -19,7 +19,7 @@ class ContourDisplay(Display):
     """
 
     def __init__(self, ds, subplot_shape=(1,), ds_name=None, **kwargs):
-        super().__init__(ds, subplot_shape, ds_name, **kwargs)
+        super().__init__(ds, subplot_shape, ds_name, secondary_y_allowed=False, **kwargs)
 
     def create_contour(
         self,
diff --git a/act/plotting/distributiondisplay.py b/act/plotting/distributiondisplay.py
index 92baa94ff4..a819bbd8aa 100644
--- a/act/plotting/distributiondisplay.py
+++ b/act/plotting/distributiondisplay.py
@@ -33,7 +33,7 @@ class DistributionDisplay(Display):
     """
 
     def __init__(self, ds, subplot_shape=(1,), ds_name=None, **kwargs):
-        super().__init__(ds, subplot_shape, ds_name, **kwargs)
+        super().__init__(ds, subplot_shape, ds_name, secondary_y_allowed=True, **kwargs)
 
     def set_xrng(self, xrng, subplot_index=(0,)):
         """
@@ -55,7 +55,7 @@ def set_xrng(self, xrng, subplot_index=(0,)):
         elif not hasattr(self, 'xrng') and len(self.axes.shape) == 1:
             self.xrng = np.zeros((self.axes.shape[0], 2), dtype='datetime64[D]')
 
-        self.axes[subplot_index].set_xlim(xrng)
+        self.axes[subplot_index][0].set_xlim(xrng)
         self.xrng[subplot_index, :] = np.array(xrng)
 
     def set_yrng(self, yrng, subplot_index=(0,)):
@@ -81,7 +81,7 @@ def set_yrng(self, yrng, subplot_index=(0,)):
         if yrng[0] == yrng[1]:
             yrng[1] = yrng[1] + 1
 
-        self.axes[subplot_index].set_ylim(yrng)
+        self.axes[subplot_index][0].set_ylim(yrng)
         self.yrng[subplot_index, :] = yrng
 
     def _get_data(self, dsname, fields):
@@ -163,13 +163,13 @@ def plot_stacked_bar_graph(
             # We will defaut the y direction to have the same # of bins as x
             sortby_bins = np.linspace(ydata.values.min(), ydata.values.max(), len(bins))
 
-        # Get the current plotting axis, add day/night background and plot data
+        # Get the current plotting axis
         if self.fig is None:
             self.fig = plt.figure()
-
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
 
         if sortby_field is not None:
             if 'units' in ydata.attrs:
@@ -189,26 +189,26 @@ def plot_stacked_bar_graph(
                     bins=[bins, sortby_bins],
                     **hist_kwargs)
             x_inds = (x_bins[:-1] + x_bins[1:]) / 2.0
-            self.axes[subplot_index].bar(
+            self.axes[subplot_index][0].bar(
                 x_inds,
                 my_hist[:, 0].flatten(),
                 label=(str(y_bins[0]) + ' to ' + str(y_bins[1])),
                 **kwargs,
             )
             for i in range(1, len(y_bins) - 1):
-                self.axes[subplot_index].bar(
+                self.axes[subplot_index][0].bar(
                     x_inds,
                     my_hist[:, i].flatten(),
                     bottom=my_hist[:, i - 1],
                     label=(str(y_bins[i]) + ' to ' + str(y_bins[i + 1])),
                     **kwargs,
                 )
-            self.axes[subplot_index].legend()
+            self.axes[subplot_index][0].legend()
         else:
             my_hist, bins = np.histogram(xdata.values.flatten(), bins=bins,
                                          density=density, **hist_kwargs)
             x_inds = (bins[:-1] + bins[1:]) / 2.0
-            self.axes[subplot_index].bar(x_inds, my_hist)
+            self.axes[subplot_index][0].bar(x_inds, my_hist)
 
         # Set Title
         if set_title is None:
@@ -220,9 +220,9 @@ def plot_stacked_bar_graph(
                     dt_utils.numpy_to_arm_date(self._ds[dsname].time.values[0]),
                 ]
             )
-        self.axes[subplot_index].set_title(set_title)
-        self.axes[subplot_index].set_ylabel('count')
-        self.axes[subplot_index].set_xlabel(xtitle)
+        self.axes[subplot_index][0].set_title(set_title)
+        self.axes[subplot_index][0].set_ylabel('count')
+        self.axes[subplot_index][0].set_xlabel(xtitle)
 
         return_dict = {}
         return_dict['plot_handle'] = self.axes[subplot_index]
@@ -306,13 +306,13 @@ def plot_size_distribution(
                 + 'length is equal to the field length!'
             )
 
-        # Get the current plotting axis, add day/night background and plot data
+        # Get the current plotting axis
         if self.fig is None:
             self.fig = plt.figure()
-
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
 
         # Set Title
         if set_title is None:
@@ -327,10 +327,10 @@ def plot_size_distribution(
             if time is not None:
                 t = pd.Timestamp(time)
                 set_title += ''.join([' at ', ':'.join([str(t.hour), str(t.minute), str(t.second)])])
-        self.axes[subplot_index].set_title(set_title)
-        self.axes[subplot_index].step(bins.values, xdata.values, **kwargs)
-        self.axes[subplot_index].set_xlabel(xtitle)
-        self.axes[subplot_index].set_ylabel(ytitle)
+        self.axes[subplot_index][0].set_title(set_title)
+        self.axes[subplot_index][0].step(bins.values, xdata.values, **kwargs)
+        self.axes[subplot_index][0].set_xlabel(xtitle)
+        self.axes[subplot_index][0].set_ylabel(ytitle)
 
         return self.axes[subplot_index]
 
@@ -412,8 +412,9 @@ def plot_stairstep_graph(
             self.fig = plt.figure()
 
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
 
         if sortby_field is not None:
             if 'units' in ydata.attrs:
@@ -433,26 +434,26 @@ def plot_stairstep_graph(
                     **hist_kwargs
                 )
             x_inds = (x_bins[:-1] + x_bins[1:]) / 2.0
-            self.axes[subplot_index].step(
+            self.axes[subplot_index][0].step(
                 x_inds,
                 my_hist[:, 0].flatten(),
                 label=(str(y_bins[0]) + ' to ' + str(y_bins[1])),
                 **kwargs,
             )
             for i in range(1, len(y_bins) - 1):
-                self.axes[subplot_index].step(
+                self.axes[subplot_index][0].step(
                     x_inds,
                     my_hist[:, i].flatten(),
                     label=(str(y_bins[i]) + ' to ' + str(y_bins[i + 1])),
                     **kwargs,
                 )
-            self.axes[subplot_index].legend()
+            self.axes[subplot_index][0].legend()
         else:
             my_hist, bins = np.histogram(xdata.values.flatten(), bins=bins,
                                          density=density, **hist_kwargs)
 
             x_inds = (bins[:-1] + bins[1:]) / 2.0
-            self.axes[subplot_index].step(x_inds, my_hist, **kwargs)
+            self.axes[subplot_index][0].step(x_inds, my_hist, **kwargs)
 
         # Set Title
         if set_title is None:
@@ -464,9 +465,9 @@ def plot_stairstep_graph(
                     dt_utils.numpy_to_arm_date(self._ds[dsname].time.values[0]),
                 ]
             )
-        self.axes[subplot_index].set_title(set_title)
-        self.axes[subplot_index].set_ylabel('count')
-        self.axes[subplot_index].set_xlabel(xtitle)
+        self.axes[subplot_index][0].set_title(set_title)
+        self.axes[subplot_index][0].set_ylabel('count')
+        self.axes[subplot_index][0].set_xlabel(xtitle)
 
         return_dict = {}
         return_dict['plot_handle'] = self.axes[subplot_index]
@@ -568,10 +569,10 @@ def plot_heatmap(
         # Get the current plotting axis, add day/night background and plot data
         if self.fig is None:
             self.fig = plt.figure()
-
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
 
         if 'units' in ydata.attrs:
             ytitle = ''.join(['(', ydata.attrs['units'], ')'])
@@ -597,7 +598,7 @@ def plot_heatmap(
         x_inds = (x_bins[:-1] + x_bins[1:]) / 2.0
         y_inds = (y_bins[:-1] + y_bins[1:]) / 2.0
         xi, yi = np.meshgrid(x_inds, y_inds, indexing='ij')
-        mesh = self.axes[subplot_index].pcolormesh(xi, yi, my_hist, shading=set_shading, **kwargs)
+        mesh = self.axes[subplot_index][0].pcolormesh(xi, yi, my_hist, shading=set_shading, **kwargs)
 
         # Set Title
         if set_title is None:
@@ -608,13 +609,13 @@ def plot_heatmap(
                     dt_utils.numpy_to_arm_date(self._ds[dsname].time.values[0]),
                 ]
             )
-        self.axes[subplot_index].set_title(set_title)
-        self.axes[subplot_index].set_ylabel(ytitle)
-        self.axes[subplot_index].set_xlabel(xtitle)
+        self.axes[subplot_index][0].set_title(set_title)
+        self.axes[subplot_index][0].set_ylabel(ytitle)
+        self.axes[subplot_index][0].set_xlabel(xtitle)
         self.add_colorbar(mesh, title='count', subplot_index=subplot_index)
 
         return_dict = {}
-        return_dict['plot_handle'] = self.axes[subplot_index]
+        return_dict['plot_handle'] = self.axes[subplot_index][0]
         return_dict['x_bins'] = x_bins
         return_dict['y_bins'] = y_bins
         return_dict['histogram'] = my_hist
@@ -634,9 +635,9 @@ def set_ratio_line(self, subplot_index=(0, )):
         if self.axes is None:
             raise RuntimeError('set_ratio_line requires the plot to be displayed.')
         # Define the xticks of the figure
-        xlims = self.axes[subplot_index].get_xticks()
+        xlims = self.axes[subplot_index][0].get_xticks()
         ratio = np.linspace(xlims[0], xlims[-1])
-        self.axes[subplot_index].plot(ratio, ratio, 'k--')
+        self.axes[subplot_index][0].plot(ratio, ratio, 'k--')
 
     def plot_scatter(self,
                      x_field,
@@ -713,15 +714,12 @@ def plot_scatter(self,
 
         # Define the axes for the figure
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
 
         # Display the scatter plot, pass keyword args for unspecified attributes
-        scc = self.axes[subplot_index].scatter(xdata,
-                                               ydata,
-                                               c=mdata,
-                                               **kwargs
-                                               )
+        scc = self.axes[subplot_index][0].scatter(xdata, ydata, c=mdata, **kwargs)
 
         # Set Title
         if set_title is None:
@@ -748,9 +746,9 @@ def plot_scatter(self,
             cbar.ax.set_ylabel(ztitle)
 
         # Define the axe title, x-axis label, y-axis label
-        self.axes[subplot_index].set_title(set_title)
-        self.axes[subplot_index].set_ylabel(ytitle)
-        self.axes[subplot_index].set_xlabel(xtitle)
+        self.axes[subplot_index][0].set_title(set_title)
+        self.axes[subplot_index][0].set_ylabel(ytitle)
+        self.axes[subplot_index][0].set_xlabel(xtitle)
 
         return self.axes[subplot_index]
 
@@ -818,8 +816,9 @@ def plot_violin(self,
 
         # Define the axes for the figure
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
 
         # Define the axe label. If units are avaiable, plot.
         if 'units' in ndata.attrs:
@@ -828,14 +827,14 @@ def plot_violin(self,
             axtitle = field
 
         # Display the scatter plot, pass keyword args for unspecified attributes
-        scc = self.axes[subplot_index].violinplot(ndata,
-                                                  positions=positions,
-                                                  vert=vert,
-                                                  showmeans=showmeans,
-                                                  showmedians=showmedians,
-                                                  showextrema=showextrema,
-                                                  **kwargs
-                                                  )
+        scc = self.axes[subplot_index][0].violinplot(ndata,
+                                                     positions=positions,
+                                                     vert=vert,
+                                                     showmeans=showmeans,
+                                                     showmedians=showmedians,
+                                                     showextrema=showextrema,
+                                                     **kwargs
+                                                     )
         if showmeans is True:
             scc['cmeans'].set_edgecolor('red')
             scc['cmeans'].set_label('mean')
@@ -853,14 +852,14 @@ def plot_violin(self,
             )
 
         # Define the axe title, x-axis label, y-axis label
-        self.axes[subplot_index].set_title(set_title)
+        self.axes[subplot_index][0].set_title(set_title)
         if vert is True:
-            self.axes[subplot_index].set_ylabel(axtitle)
+            self.axes[subplot_index][0].set_ylabel(axtitle)
             if positions is None:
-                self.axes[subplot_index].set_xticks([])
+                self.axes[subplot_index][0].set_xticks([])
         else:
-            self.axes[subplot_index].set_xlabel(axtitle)
+            self.axes[subplot_index][0].set_xlabel(axtitle)
             if positions is None:
-                self.axes[subplot_index].set_yticks([])
+                self.axes[subplot_index][0].set_yticks([])
 
-        return self.axes[subplot_index]
+        return self.axes[subplot_index][0]
diff --git a/act/plotting/geodisplay.py b/act/plotting/geodisplay.py
index f26379df39..1356aef231 100644
--- a/act/plotting/geodisplay.py
+++ b/act/plotting/geodisplay.py
@@ -44,7 +44,7 @@ def __init__(self, ds, ds_name=None, **kwargs):
             raise ImportError(
                 'Cartopy needs to be installed on your ' 'system to make geographic display plots.'
             )
-        super().__init__(ds, ds_name, **kwargs)
+        super().__init__(ds, ds_name, secondary_y_allowed=False, **kwargs)
         if self.fig is None:
             self.fig = plt.figure(**kwargs)
 
diff --git a/act/plotting/plot.py b/act/plotting/plot.py
index 76ad2693e7..2dccff414b 100644
--- a/act/plotting/plot.py
+++ b/act/plotting/plot.py
@@ -66,12 +66,16 @@ class with this set to None will create a new figure handle. See the
         loaded.
     subplot_kw : dict, optional
         The kwargs to pass into :func:`fig.subplots`
+    secondary_y_allowed : boolean
+        If the plot type allows a secondary y axis
+
     **kwargs : keywords arguments
         Keyword arguments passed to :func:`plt.figure`.
 
     """
 
-    def __init__(self, ds, subplot_shape=(1,), ds_name=None, subplot_kw=None, **kwargs):
+    def __init__(self, ds, subplot_shape=(1,), ds_name=None, subplot_kw=None,
+                 secondary_y_allowed=True, **kwargs):
         if isinstance(ds, xr.Dataset):
             if 'datastream' in ds.attrs.keys() is not None:
                 self._ds = {ds.attrs['datastream']: ds}
@@ -119,9 +123,11 @@ def __init__(self, ds, subplot_shape=(1,), ds_name=None, subplot_kw=None, **kwar
         self.plot_vars = []
         self.cbs = []
         if subplot_shape is not None:
-            self.add_subplots(subplot_shape, subplot_kw=subplot_kw, **kwargs)
+            self.add_subplots(subplot_shape, subplot_kw=subplot_kw,
+                              secondary_y_allowed=secondary_y_allowed, **kwargs)
 
-    def add_subplots(self, subplot_shape=(1,), subplot_kw=None, **kwargs):
+    def add_subplots(self, subplot_shape=(1,), secondary_y=False, subplot_kw=None,
+                     secondary_y_allowed=True, **kwargs):
         """
         Adds subplots to the Display object. The current
         figure in the object will be deleted and overwritten.
@@ -159,10 +165,29 @@ def add_subplots(self, subplot_shape=(1,), subplot_kw=None, **kwargs):
             self.yrng = np.zeros((subplot_shape[0], 2))
         else:
             raise ValueError('subplot_shape must be a 1 or 2 dimensional' + 'tuple list, or array!')
+        # Create dummy ax to add secondary y-axes to
+        if secondary_y_allowed:
+            dummy_ax = np.empty(list(ax.shape) + [2], dtype=plt.Axes)
+            for i, axis in enumerate(dummy_ax):
+                if len(axis.shape) == 1:
+                    dummy_ax[i, 0] = ax[i]
+                    try:
+                        dummy_ax[i, 1] = ax[i].twinx()
+                    except Exception:
+                        dummy_ax[i, 1] = None
+                else:
+                    for j, axis2 in enumerate(axis):
+                        dummy_ax[i, j, 0] = ax[i, j]
+                        try:
+                            dummy_ax[i, j, 1] = ax[i, j].twinx()
+                        except Exception:
+                            dummy_ax[i, j, 1] = None
+        else:
+            dummy_ax = ax
         self.fig = fig
-        self.axes = ax
+        self.axes = dummy_ax
 
-    def put_display_in_subplot(self, display, subplot_index):
+    def put_display_in_subplot(self, display, subplot_index, y_axis_index=0):
         """
         This will place a Display object into a specific subplot.
         The display object must only have one subplot.
@@ -186,11 +211,17 @@ def put_display_in_subplot(self, display, subplot_index):
             raise RuntimeError(
                 'Only single plots can be made as subplots ' + 'of another Display object!'
             )
+        if len(np.shape(display.axes)) == 1:
+            my_projection = display.axes[0].name
+        else:
+            my_projection = display.axes[0][y_axis_index].name
 
-        my_projection = display.axes[0].name
         plt.close(display.fig)
         display.fig = self.fig
-        self.fig.delaxes(self.axes[subplot_index])
+        if len(np.shape(self.axes)) == 1:
+            self.fig.delaxes(self.axes[subplot_index])
+        else:
+            self.fig.delaxes(self.axes[subplot_index][y_axis_index])
         the_shape = self.axes.shape
         if len(the_shape) == 1:
             second_value = 1
@@ -261,7 +292,10 @@ def add_colorbar(self, mappable, title=None, subplot_index=(0,), pad=None,
             raise RuntimeError('add_colorbar requires the plot ' 'to be displayed.')
 
         fig = self.fig
-        ax = self.axes[subplot_index]
+        if np.size(self.axes[subplot_index]) > 1:
+            ax = self.axes[subplot_index][0]
+        else:
+            ax = self.axes[subplot_index]
 
         if pad is None:
             pad = 0.01
@@ -347,6 +381,9 @@ def plot_group(self, func_name, dsname=None, **kwargs):
             raise RuntimeError("The specified string is not a function of "
                                "the Display object.")
         subplot_shape = self.display.axes.shape
+        if len(subplot_shape) > 2:
+            subplot_shape = subplot_shape[0:-1]
+
         i = 0
         wrap_around = False
         old_ds = self.display._ds
@@ -362,7 +399,7 @@ def plot_group(self, func_name, dsname=None, **kwargs):
                     if len(subplot_shape) == 2:
                         subplot_index = (int(i / subplot_shape[1]), i % subplot_shape[1])
                     else:
-                        subplot_index = (i % subplot_shape[0],)
+                        subplot_index = (i % subplot_shape[0], 0)
                     args, varargs, varkw, _, _, _, _ = inspect.getfullargspec(func)
                     if "subplot_index" in args:
                         kwargs["subplot_index"] = subplot_index
@@ -396,25 +433,33 @@ def plot_group(self, func_name, dsname=None, **kwargs):
                 if len(subplot_shape) == 2:
                     subplot_index = (int(i / subplot_shape[1]), i % subplot_shape[1])
                 else:
-                    subplot_index = (i % subplot_shape[0],)
-                self.display.axes[subplot_index].axis('off')
+                    subplot_index = (i % subplot_shape[0], 0)
+                if np.size(self.display.axes) == 1:
+                    self.display.axes[subplot_index].axis('off')
+                elif np.size(self.display.axes[subplot_index]) > 1:
+                    self.display.axes[subplot_index][0].axis('off')
+                    self.display.axes[subplot_index][1].axis('off')
+                else:
+                    self.display.axes[subplot_index].axis('off')
                 i = i + 1
 
         for i in range(1, np.prod(subplot_shape)):
             if len(subplot_shape) == 2:
                 subplot_index = (int(i / subplot_shape[1]), i % subplot_shape[1])
             else:
-                subplot_index = (i % subplot_shape[0],)
+                subplot_index = (i % subplot_shape[0], 0)
             try:
                 self.display.axes[subplot_index].get_legend().remove()
             except AttributeError:
                 pass
-
         if self.isTimeSeriesDisplay:
             key_list = list(self.display._ds.keys())
             for k in key_list:
                 time_min, time_max = self.xlims[k]
-                subplot_index = self.mapping[k]
+                if len(self.mapping[k]) == 1:
+                    subplot_index = self.mapping[k] + (0,)
+                else:
+                    subplot_index = self.mapping[k]
                 self.display.set_xrng([time_min, time_max], subplot_index)
 
         self.display._ds = old_ds
diff --git a/act/plotting/skewtdisplay.py b/act/plotting/skewtdisplay.py
index 2203d3cf9c..c70875dc86 100644
--- a/act/plotting/skewtdisplay.py
+++ b/act/plotting/skewtdisplay.py
@@ -56,7 +56,8 @@ def __init__(self, ds, subplot_shape=(1,), subplot=None, ds_name=None, set_fig=N
         # We want to use our routine to handle subplot adding, not the main
         # one
         new_kwargs = kwargs.copy()
-        super().__init__(ds, None, ds_name, subplot_kw=dict(projection='skewx'), **new_kwargs)
+        super().__init__(ds, None, ds_name, subplot_kw=dict(projection='skewx'),
+                         secondary_y_allowed=False, **new_kwargs)
 
         # Make a SkewT object for each subplot
         self.add_subplots(subplot_shape, set_fig=set_fig, subplot=subplot, **kwargs)
diff --git a/act/plotting/timeseriesdisplay.py b/act/plotting/timeseriesdisplay.py
index 948c48028d..4ecc43c80d 100644
--- a/act/plotting/timeseriesdisplay.py
+++ b/act/plotting/timeseriesdisplay.py
@@ -51,8 +51,8 @@ class TimeSeriesDisplay(Display):
 
     """
 
-    def __init__(self, ds, subplot_shape=(1,), ds_name=None, **kwargs):
-        super().__init__(ds, subplot_shape, ds_name, **kwargs)
+    def __init__(self, ds, subplot_shape=(1,), ds_name=None, secondary_y_allowed=True, **kwargs):
+        super().__init__(ds, subplot_shape, ds_name, secondary_y_allowed=True, **kwargs)
 
     def day_night_background(self, dsname=None, subplot_index=(0,)):
         """
@@ -94,7 +94,11 @@ def day_night_background(self, dsname=None, subplot_index=(0,)):
         if self.axes is None:
             raise RuntimeError('day_night_background requires the plot to ' 'be displayed.')
 
-        ax = self.axes[subplot_index]
+        # Default to the left axis
+        if np.size(self.axes[subplot_index]) > 1:
+            ax = self.axes[subplot_index][0]
+        else:
+            ax = self.axes[subplot_index]
 
         # Find variable names for latitude and longitude
         variables = list(self._ds[dsname].data_vars)
@@ -196,7 +200,7 @@ def day_night_background(self, dsname=None, subplot_index=(0,)):
         for ii in noon:
             ax.axvline(x=ii, linestyle='--', color='y', zorder=1)
 
-    def set_xrng(self, xrng, subplot_index=(0,)):
+    def set_xrng(self, xrng, subplot_index=(0, 0), y_axis_index=0):
         """
         Sets the x range of the plot.
 
@@ -224,8 +228,10 @@ def set_xrng(self, xrng, subplot_index=(0,)):
                       'Expanding range by 2 seconds.\n')
                 xrng[0] -= dt.timedelta(seconds=1)
                 xrng[1] += dt.timedelta(seconds=1)
-
-        self.axes[subplot_index].set_xlim(xrng)
+        if np.size(self.axes[subplot_index]) > 1:
+            self.axes[subplot_index][y_axis_index].set_xlim(xrng)
+        else:
+            self.axes[subplot_index].set_xlim(xrng)
 
         # Make sure that the xrng value is a numpy array not pandas
         if isinstance(xrng[0], pd.Timestamp):
@@ -242,7 +248,7 @@ def set_xrng(self, xrng, subplot_index=(0,)):
             self.xrng[subplot_index][0] = xrng[0].astype('datetime64[D]').astype(float)
             self.xrng[subplot_index][1] = xrng[1].astype('datetime64[D]').astype(float)
 
-    def set_yrng(self, yrng, subplot_index=(0,), match_axes_ylimits=False):
+    def set_yrng(self, yrng, subplot_index=(0,), match_axes_ylimits=False, y_axis_index=0):
         """
         Sets the y range of the plot.
 
@@ -257,6 +263,8 @@ def set_yrng(self, yrng, subplot_index=(0,), match_axes_ylimits=False):
             If True, all axes in the display object will have matching
             provided ylims. Default is False. This is especially useful
             when utilizing a groupby display with many axes.
+        y_axis_index : int
+            0 = left y axis, 1 = right y axis
 
         """
         if self.axes is None:
@@ -274,9 +282,12 @@ def set_yrng(self, yrng, subplot_index=(0,), match_axes_ylimits=False):
         if match_axes_ylimits:
             for i in range(self.axes.shape[0]):
                 for j in range(self.axes.shape[1]):
-                    self.axes[i, j].set_ylim(yrng)
+                    self.axes[i, j, y_axis_index].set_ylim(yrng)
         else:
-            self.axes[subplot_index].set_ylim(yrng)
+            if np.size(self.axes[subplot_index]) > 1:
+                self.axes[subplot_index][y_axis_index].set_ylim(yrng)
+            else:
+                self.axes[subplot_index].set_ylim(yrng)
 
         try:
             self.yrng[subplot_index, :] = yrng
@@ -287,7 +298,7 @@ def plot(
         self,
         field,
         dsname=None,
-        subplot_index=(0,),
+        subplot_index=(0, ),
         cmap=None,
         set_title=None,
         add_nan=False,
@@ -393,6 +404,7 @@ def plot(
             move to right negative values move to left.
         secondary_y : boolean
             Option to plot on secondary y axis.
+            This will automatically change the color of the axis to match the line
         y_axis_flag_meanings : boolean or int
             When set to True and plotting state variable with flag_values and
             flag_meanings attributes will replace y axis numerical values
@@ -486,19 +498,27 @@ def plot(
         else:
             ydata = None
 
-        # Get the current plotting axis, add day/night background and plot data
+        # Get the current plotting axis
         if self.fig is None:
             self.fig = plt.figure()
 
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
 
-        # Set up secondary y axis if requested
         if secondary_y is False:
-            ax = self.axes[subplot_index]
+            y_axis_index = 0
+            if np.size(self.axes[subplot_index]) > 1:
+                ax = self.axes[subplot_index][y_axis_index]
+                self.axes[subplot_index][1].get_yaxis().set_visible(False)
+            else:
+                ax = self.axes[subplot_index]
         else:
-            ax = self.axes[subplot_index].twinx()
+            y_axis_index = 1
+            ax = self.axes[subplot_index][y_axis_index]
+            self.axes[subplot_index][1].get_yaxis().set_visible(True)
+            match_line_label_color = True
 
         if colorbar_labels is not None:
             flag_values = list(colorbar_labels.keys())
@@ -656,6 +676,11 @@ def plot(
         if not y_axis_flag_meanings:
             if match_line_label_color and len(ax.get_lines()) > 0:
                 ax.set_ylabel(ytitle, color=ax.get_lines()[0].get_color())
+                ax.tick_params(axis='y', colors=ax.get_lines()[0].get_color())
+                if y_axis_index == 0:
+                    ax.spines['left'].set_color(ax.get_lines()[0].get_color())
+                if y_axis_index == 1:
+                    ax.spines['right'].set_color(ax.get_lines()[0].get_color())
             else:
                 ax.set_ylabel(ytitle)
 
@@ -715,12 +740,7 @@ def plot(
                 if yrng[1] > current_yrng[1]:
                     yrng[1] = current_yrng[1]
 
-            # Set y range the normal way if not secondary y
-            # If secondary, just use set_ylim
-            if secondary_y is False:
-                self.set_yrng(yrng, subplot_index)
-            else:
-                ax.set_ylim(yrng)
+            self.set_yrng(yrng, subplot_index, y_axis_index=y_axis_index)
 
         # Set X Format
         if len(subplot_index) == 1:
@@ -762,7 +782,6 @@ def plot(
                 self.add_colorbar(
                     mesh, title=cbar_title, subplot_index=subplot_index, pad=cbar_h_adjust
                 )
-
         return ax
 
     def plot_barbs_from_spd_dir(
@@ -969,16 +988,26 @@ def plot_barbs_from_u_v(
         if self.fig is None:
             self.fig = plt.figure()
 
+        # Set up or get current axes
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
+
+        # Setting up in case there is a use case in the future for a secondary y
+        y_axis_index = 0
+
+        if len(np.shape(self.axes)) == 1:
+            ax = self.axes[subplot_index]
+        else:
+            ax = self.axes[subplot_index][y_axis_index]
 
         if ydata is None:
             ydata = np.ones(xdata.shape)
             if 'cmap' in kwargs.keys():
                 map_color = np.sqrt(np.power(u[::barb_step_x], 2) + np.power(v[::barb_step_x], 2))
                 map_color[np.isnan(map_color)] = 0
-                ax = self.axes[subplot_index].barbs(
+                barbs = ax.barbs(
                     xdata[::barb_step_x],
                     ydata[::barb_step_x],
                     u[::barb_step_x],
@@ -987,20 +1016,20 @@ def plot_barbs_from_u_v(
                     **kwargs,
                 )
                 plt.colorbar(
-                    ax,
-                    ax=[self.axes[subplot_index]],
+                    barbs,
+                    ax=[ax],
                     label='Wind Speed (' + self._ds[dsname][u_field].attrs['units'] + ')',
                 )
 
             else:
-                self.axes[subplot_index].barbs(
+                ax.barbs(
                     xdata[::barb_step_x],
                     ydata[::barb_step_x],
                     u[::barb_step_x],
                     v[::barb_step_x],
                     **kwargs,
                 )
-            self.axes[subplot_index].set_yticks([])
+            ax.set_yticks([])
 
         else:
             if 'cmap' in kwargs.keys():
@@ -1009,7 +1038,7 @@ def plot_barbs_from_u_v(
                     + np.power(v[::barb_step_x, ::barb_step_y], 2)
                 )
                 map_color[np.isnan(map_color)] = 0
-                ax = self.axes[subplot_index].barbs(
+                barbs = ax.barbs(
                     xdata[::barb_step_x, ::barb_step_y],
                     ydata[::barb_step_x, ::barb_step_y],
                     u[::barb_step_x, ::barb_step_y],
@@ -1018,12 +1047,12 @@ def plot_barbs_from_u_v(
                     **kwargs,
                 )
                 plt.colorbar(
-                    ax,
-                    ax=[self.axes[subplot_index]],
+                    barbs,
+                    ax=[ax],
                     label='Wind Speed (' + self._ds[dsname][u_field].attrs['units'] + ')',
                 )
             else:
-                ax = self.axes[subplot_index].barbs(
+                barbs = ax.barbs(
                     xdata[::barb_step_x, ::barb_step_y],
                     ydata[::barb_step_x, ::barb_step_y],
                     u[::barb_step_x, ::barb_step_y],
@@ -1044,11 +1073,11 @@ def plot_barbs_from_u_v(
                 ]
             )
 
-        self.axes[subplot_index].set_title(set_title)
+        ax.set_title(set_title)
 
         # Set YTitle
         if 'ytitle' in locals():
-            self.axes[subplot_index].set_ylabel(ytitle)
+            ax.set_ylabel(ytitle)
 
         # Set X Limit - We want the same time axes for all subplots
         time_rng = [xdata.min(), xdata.max()]
@@ -1085,10 +1114,14 @@ def plot_barbs_from_u_v(
 
         # Put on an xlabel, but only if we are making the bottom-most plot
         if subplot_index[0] == self.axes.shape[0] - 1:
-            self.axes[subplot_index].set_xlabel('Time [UTC]')
+            ax.set_xlabel('Time [UTC]')
 
         myFmt = common.get_date_format(days)
-        self.axes[subplot_index].xaxis.set_major_formatter(myFmt)
+        ax.xaxis.set_major_formatter(myFmt)
+        if len(np.shape(self.axes)) == 1:
+            self.axes[subplot_index] = ax
+        else:
+            self.axes[subplot_index][y_axis_index] = ax
         return self.axes[subplot_index]
 
     def plot_time_height_xsection_from_1d_data(
@@ -1186,11 +1219,21 @@ def plot_time_height_xsection_from_1d_data(
         if self.fig is None:
             self.fig = plt.figure()
 
+        # Set up or get current axes
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
+
+        # Setting up in case there is a use case in the future for a secondary y
+        y_axis_index = 0
 
-        mesh = self.axes[subplot_index].pcolormesh(
+        if len(np.shape(self.axes)) == 1:
+            ax = self.axes[subplot_index]
+        else:
+            ax = self.axes[subplot_index][y_axis_index]
+
+        mesh = ax.pcolormesh(
             x_times, y_levels, np.transpose(data), shading=set_shading, **kwargs
         )
 
@@ -1207,11 +1250,11 @@ def plot_time_height_xsection_from_1d_data(
                 ]
             )
 
-        self.axes[subplot_index].set_title(set_title)
+        ax.set_title(set_title)
 
         # Set YTitle
         if 'ytitle' in locals():
-            self.axes[subplot_index].set_ylabel(ytitle)
+            ax.set_ylabel(ytitle)
 
         # Set X Limit - We want the same time axes for all subplots
         time_rng = [x_times[0], x_times[-1]]
@@ -1248,7 +1291,7 @@ def plot_time_height_xsection_from_1d_data(
 
         # Put on an xlabel, but only if we are making the bottom-most plot
         if subplot_index[0] == self.axes.shape[0] - 1:
-            self.axes[subplot_index].set_xlabel('Time [UTC]')
+            ax.set_xlabel('Time [UTC]')
 
         if ydata is not None:
             if cbar_label is None:
@@ -1256,7 +1299,7 @@ def plot_time_height_xsection_from_1d_data(
             else:
                 self.add_colorbar(mesh, title=cbar_label, subplot_index=subplot_index)
         myFmt = common.get_date_format(days)
-        self.axes[subplot_index].xaxis.set_major_formatter(myFmt)
+        ax.xaxis.set_major_formatter(myFmt)
 
         return self.axes[subplot_index]
 
@@ -1442,10 +1485,15 @@ def qc_flag_block_plot(
 
         # Set up or get current axes
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
+
+        # Setting y_axis_index in case there is a use case in the future
+        # to plot the QC on the secondary y-axis
+        y_axis_index = 0
 
-        ax = self.axes[subplot_index]
+        ax = self.axes[subplot_index][y_axis_index]
 
         # Set X Limit - We want the same time axes for all subplots
         data = self._ds[dsname][data_field]
@@ -1719,14 +1767,17 @@ def fill_between(
             self.fig = plt.figure()
 
         if self.axes is None:
-            self.axes = np.array([plt.axes()])
-            self.fig.add_axes(self.axes[0])
+            self.axes = np.array([[plt.axes(), plt.axes().twinx()]])
+            for a in self.axes[0]:
+                self.fig.add_axes(a)
 
         # Set ax to appropriate axis
         if secondary_y is False:
-            ax = self.axes[subplot_index]
+            y_axis_index = 0
+            ax = self.axes[subplot_index][y_axis_index]
         else:
-            ax = self.axes[subplot_index].twinx()
+            y_axis_index = 1
+            ax = self.axes[subplot_index][y_axis_index]
 
         ax.fill_between(xdata.values, data, **kwargs)
 
@@ -1748,7 +1799,7 @@ def fill_between(
 
         # Put on an xlabel, but only if we are making the bottom-most plot
         if subplot_index[0] == self.axes.shape[0] - 1:
-            self.axes[subplot_index].set_xlabel('Time [UTC]')
+            ax.set_xlabel('Time [UTC]')
 
         # Set YTitle
         ax.set_ylabel(ytitle)
@@ -1765,5 +1816,5 @@ def fill_between(
             )
         if secondary_y is False:
             ax.set_title(set_title)
-
+        self.axes[subplot_index][y_axis_index] = ax
         return self.axes[subplot_index]
diff --git a/act/plotting/windrosedisplay.py b/act/plotting/windrosedisplay.py
index 365cc1325f..4d7b561e33 100644
--- a/act/plotting/windrosedisplay.py
+++ b/act/plotting/windrosedisplay.py
@@ -36,7 +36,8 @@ class and has therefore has the same attributes as that class.
     """
 
     def __init__(self, ds, subplot_shape=(1,), ds_name=None, **kwargs):
-        super().__init__(ds, subplot_shape, ds_name, subplot_kw=dict(projection='polar'), **kwargs)
+        super().__init__(ds, subplot_shape, ds_name, subplot_kw=dict(projection='polar'),
+                         secondary_y_allowed=False, **kwargs)
 
     def set_thetarng(self, trng=(0.0, 360.0), subplot_index=(0,)):
         """
@@ -194,8 +195,18 @@ def plot(
         our_cmap = matplotlib.colormaps.get_cmap(cmap)
         our_colors = our_cmap(np.linspace(0, 1, len(spd_bins)))
 
+        # Make sure we're dealing with the right axes style
+        if len(np.shape(self.axes)) == 1:
+            ax = self.axes[subplot_index]
+            if np.size(ax) > 1:
+                ax = ax[0]
+        elif len(self.axes[subplot_index]) == 2:
+            ax = self.axes[subplot_index][0]
+        else:
+            ax = self.axes[subplot_index]
+
         bars = [
-            self.axes[subplot_index].bar(
+            ax.bar(
                 mins,
                 wind_hist[:, 0],
                 bottom=0,
@@ -210,7 +221,7 @@ def plot(
             # Changing the bottom to be a sum of the previous speeds so that
             # it positions it correctly - Adam Theisen
             bars.append(
-                self.axes[subplot_index].bar(
+                ax.bar(
                     mins,
                     wind_hist[:, i],
                     label=the_label,
@@ -220,16 +231,16 @@ def plot(
                     **kwargs,
                 )
             )
-        self.axes[subplot_index].legend(
+        ax.legend(
             loc=legend_loc, bbox_to_anchor=legend_bbox, title=legend_title
         )
-        self.axes[subplot_index].set_theta_zero_location('N')
-        self.axes[subplot_index].set_theta_direction(-1)
+        ax.set_theta_zero_location('N')
+        ax.set_theta_direction(-1)
 
         # Add an annulus with text stating % of time calm
         pct_calm = np.sum(spd_data <= calm_threshold) / len(spd_data) * 100
-        self.axes[subplot_index].set_rorigin(-2.5)
-        self.axes[subplot_index].annotate(
+        ax.set_rorigin(-2.5)
+        ax.annotate(
             '%3.2f%%\n calm' % pct_calm, xy=(0, -2.5), ha='center', va='center'
         )
 
@@ -237,8 +248,8 @@ def plot(
         tick_max = tick_interval * round(np.nanmax(np.cumsum(wind_hist, axis=1)) / tick_interval)
         rticks = np.arange(0, tick_max, tick_interval)
         rticklabels = [('%d' % x + '%') for x in rticks]
-        self.axes[subplot_index].set_rticks(rticks)
-        self.axes[subplot_index].set_yticklabels(rticklabels)
+        ax.set_rticks(rticks)
+        ax.set_yticklabels(rticklabels)
 
         # Set Title
         if set_title is None:
@@ -249,8 +260,16 @@ def plot(
                     dt_utils.numpy_to_arm_date(self._ds[dsname].time.values[0]),
                 ]
             )
-        self.axes[subplot_index].set_title(set_title)
-        return self.axes[subplot_index]
+        ax.set_title(set_title)
+
+        if len(np.shape(self.axes)) == 1:
+            self.axes[subplot_index] = ax
+        elif len(self.axes[subplot_index]) == 2:
+            self.axes[subplot_index][0] = ax
+        else:
+            self.axes[subplot_index] = ax
+
+        return ax
 
     def plot_data(
         self,
diff --git a/act/plotting/xsectiondisplay.py b/act/plotting/xsectiondisplay.py
index d989594c0c..86e2a8ef51 100644
--- a/act/plotting/xsectiondisplay.py
+++ b/act/plotting/xsectiondisplay.py
@@ -72,7 +72,7 @@ class and has therefore has the same attributes as that class.
     """
 
     def __init__(self, ds, subplot_shape=(1,), ds_name=None, **kwargs):
-        super().__init__(ds, subplot_shape, ds_name, **kwargs)
+        super().__init__(ds, subplot_shape, ds_name, secondary_y_allowed=False, **kwargs)
 
     def set_subplot_to_map(self, subplot_index):
         total_num_plots = self.axes.shape
diff --git a/act/tests/baseline/test_secondary_y.png b/act/tests/baseline/test_secondary_y.png
new file mode 100644
index 0000000000..4ce523778a
Binary files /dev/null and b/act/tests/baseline/test_secondary_y.png differ
diff --git a/act/tests/test_plotting.py b/act/tests/test_plotting.py
index b549594137..8c4075ab6d 100644
--- a/act/tests/test_plotting.py
+++ b/act/tests/test_plotting.py
@@ -1,6 +1,5 @@
 import glob
 from datetime import datetime
-
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
@@ -55,9 +54,10 @@ def test_plot():
     windrose = WindRoseDisplay(met)
     display.put_display_in_subplot(windrose, subplot_index=(1, 1))
     windrose.plot('wdir_vec_mean', 'wspd_vec_mean', spd_bins=np.linspace(0, 10, 4))
-    windrose.axes[0].legend(loc='best')
+    windrose.axes[0, 0].legend(loc='best')
     met.close()
 
+    return display.fig
     try:
         return display.fig
     finally:
@@ -389,6 +389,19 @@ def test_xsection_plot():
 @pytest.mark.mpl_image_compare(tolerance=30)
 def test_xsection_plot_map():
     radar_ds = arm.read_netcdf(sample_files.EXAMPLE_VISST, combine='nested', concat_dim='time')
+    xsection = XSectionDisplay(radar_ds, figsize=(15, 8))
+    xsection.plot_xsection_map(
+        None,
+        'ir_temperature',
+        vmin=220,
+        vmax=300,
+        cmap='Greys',
+        x='longitude',
+        y='latitude',
+        isel_kwargs={'time': 0},
+    )
+    radar_ds.close()
+    return xsection.fig
 
     try:
         xsection = XSectionDisplay(radar_ds, figsize=(15, 8))
@@ -1202,7 +1215,7 @@ def test_match_ylimits_plot():
     display = act.plotting.TimeSeriesDisplay(ds, figsize=(10, 8), subplot_shape=(2, 2))
     groupby = display.group_by('day')
     groupby.plot_group('plot', None, field='temp_mean', marker=' ')
-    groupby.display.set_yrng([0, 20], match_axes_ylimits=True)
+    groupby.display.set_yrng([-20, 20], match_axes_ylimits=True)
     ds.close()
     return display.fig
 
@@ -1323,3 +1336,15 @@ def test_scatter():
     ds.close()
 
     return display.fig
+
+
+@pytest.mark.mpl_image_compare(tolerance=30)
+def test_secondary_y():
+    ds = act.io.armfiles.read_netcdf(sample_files.EXAMPLE_MET1)
+    display = act.plotting.TimeSeriesDisplay(ds, figsize=(10, 6))
+    display.plot('temp_mean', match_line_label_color=True)
+    display.plot('rh_mean', secondary_y=True, color='orange')
+    display.day_night_background()
+    ds.close()
+
+    return display.fig
diff --git a/examples/discovery/plot_asos_temp.py b/examples/discovery/plot_asos_temp.py
index ac228576ff..ef1235e4b7 100644
--- a/examples/discovery/plot_asos_temp.py
+++ b/examples/discovery/plot_asos_temp.py
@@ -19,5 +19,5 @@
 display = act.plotting.TimeSeriesDisplay(my_asoses['ORD'], subplot_shape=(2,), figsize=(15, 10))
 display.plot('temp', subplot_index=(0,))
 display.plot_barbs_from_u_v(u_field='u', v_field='v', subplot_index=(1,))
-display.axes[1].set_ylim([0, 2])
+display.axes[1, 0].set_ylim([0, 2])
 plt.show()
diff --git a/examples/discovery/plot_noaa_fmcw_moment.py b/examples/discovery/plot_noaa_fmcw_moment.py
index 50bad32ead..25f3d1c05b 100644
--- a/examples/discovery/plot_noaa_fmcw_moment.py
+++ b/examples/discovery/plot_noaa_fmcw_moment.py
@@ -55,5 +55,5 @@
     subplot_index=(1,),
 )
 # Adjust ylims of parsivel plot.
-display.axes[1].set_ylim([0, 10])
+display.axes[1, 0].set_ylim([0, 10])
 plt.show()
diff --git a/examples/plotting/plot_hist_kwargs.py b/examples/plotting/plot_hist_kwargs.py
index 481198076a..477d4414b4 100644
--- a/examples/plotting/plot_hist_kwargs.py
+++ b/examples/plotting/plot_hist_kwargs.py
@@ -18,7 +18,7 @@
 
 # Plot data
 hist_kwargs = {'range': (-10, 10)}
-histdisplay = act.plotting.HistogramDisplay(met_ds)
+histdisplay = act.plotting.DistributionDisplay(met_ds)
 histdisplay.plot_stacked_bar_graph('temp_mean', bins=np.arange(-40, 40, 5),
                                    hist_kwargs=hist_kwargs)
 plt.show()
diff --git a/examples/plotting/plot_scatter.py b/examples/plotting/plot_scatter.py
index c221e49728..97c7faafb6 100644
--- a/examples/plotting/plot_scatter.py
+++ b/examples/plotting/plot_scatter.py
@@ -43,19 +43,19 @@
 p = np.poly1d(z)
 
 # Plot the best fit line
-display.axes[0].plot(ds['true_airspeed'],
-                     p(ds['true_airspeed']),
-                     'r',
-                     linewidth=2
-                     )
+display.axes[0, 0].plot(ds['true_airspeed'],
+                        p(ds['true_airspeed']),
+                        'r',
+                        linewidth=2
+                        )
 
 # Display the line equation
-display.axes[0].text(45,
-                     135,
-                     "y = %.3fx + (%.3f)" % (z[0], z[1]),
-                     color='r',
-                     fontsize=12
-                     )
+display.axes[0, 0].text(45,
+                        135,
+                        "y = %.3fx + (%.3f)" % (z[0], z[1]),
+                        color='r',
+                        fontsize=12
+                        )
 
 # Calculate Pearson Correlation Coefficient
 cc_conc = pearsonr(ds['true_airspeed'],
@@ -63,18 +63,18 @@
                    )
 
 # Display the Pearson CC
-display.axes[0].text(45,
-                     130,
-                     "Pearson CC: %.2f" % (cc_conc[0]),
-                     fontsize=12
-                     )
+display.axes[0, 0].text(45,
+                        130,
+                        "Pearson CC: %.2f" % (cc_conc[0]),
+                        fontsize=12
+                        )
 
 # Display the total number of samples
-display.axes[0].text(45,
-                     125,
-                     "N = %.0f" % (ds['true_airspeed'].data.shape[0]),
-                     fontsize=12
-                     )
+display.axes[0, 0].text(45,
+                        125,
+                        "N = %.0f" % (ds['true_airspeed'].data.shape[0]),
+                        fontsize=12
+                        )
 
 # Display the 1:1 ratio line
 display.set_ratio_line()
diff --git a/examples/plotting/plot_secondary_y.py b/examples/plotting/plot_secondary_y.py
new file mode 100644
index 0000000000..c0e8070a16
--- /dev/null
+++ b/examples/plotting/plot_secondary_y.py
@@ -0,0 +1,49 @@
+"""
+Secondary Y-Axis Plotting
+-------------------------
+
+This example shows how to use the new capability
+to plot on the secondary y-axis. Previous versions
+of ACT only returned one axis object per plot, even
+when there was a secondary y-axis.  The new functionality
+will return two axis objects per plot for the left and
+right y axes.
+
+"""
+
+
+import act
+import matplotlib.pyplot as plt
+import xarray as xr
+
+# Read in the data from a MET file
+ds = act.io.armfiles.read_netcdf(act.tests.EXAMPLE_MET1)
+
+# Plot temperature and relative humidity with RH on the right axis
+display = act.plotting.TimeSeriesDisplay(ds, figsize=(10, 6))
+
+# Plot the data and make the y-axes color match the lines
+# Note, you need to specify the color for the secondary-y plot
+# if you want it different from the primary y-axis
+display.plot('temp_mean', match_line_label_color=True)
+display.plot('rh_mean', secondary_y=True, color='orange')
+display.day_night_background()
+
+# In a slight change, the axes returned as part of the display object
+# for TimeSeries and DistributionDisplay now return a 2D array instead
+# of a 1D array.   The second dimension is the axes object for the right
+# axis which is automatically created.
+# It can still be used like before for modifications after ACT plotting
+
+# The left axis will have an index of 0
+#               \/
+display.axes[0, 0].set_yticks([-5, 0, 5])
+display.axes[0, 0].set_yticklabels(["That's cold", "Freezing", "Above Freezing"])
+
+# The right axis will have an index of 1
+#               \/
+display.axes[0, 1].set_yticks([65, 75, 85])
+display.axes[0, 1].set_yticklabels(['Not as humid', 'Slightly Humid', 'Humid'])
+
+plt.tight_layout()
+plt.show()
diff --git a/examples/plotting/plot_violin.py b/examples/plotting/plot_violin.py
index 0ae527314d..4d4395bd5c 100644
--- a/examples/plotting/plot_violin.py
+++ b/examples/plotting/plot_violin.py
@@ -33,12 +33,9 @@
                     )
 
 # Update the tick information
-display.axes[0].set_xticks([0.5, 1, 2, 2.5])
-display.axes[0].set_xticklabels(['',
-                                 'Ambient Air\nTemp',
-                                 'Total\nTemperature',
-                                 '']
-                                )
+display.axes[0, 0].set_xticks([0.5, 1, 2, 2.5])
+ticks = ['', 'Ambient Air\nTemp', 'Total\nTemperature', '']
+display.axes[0, 0].set_xticklabels(ticks)
 
 # Update the y-axis label
-display.axes[0].set_ylabel('Temperature Observations [C]')
+display.axes[0, 0].set_ylabel('Temperature Observations [C]')
diff --git a/examples/workflows/plot_weighted_average.py b/examples/workflows/plot_weighted_average.py
index 47fd4eae90..77ae964926 100644
--- a/examples/workflows/plot_weighted_average.py
+++ b/examples/workflows/plot_weighted_average.py
@@ -90,5 +90,5 @@
 )
 display.plot('weighted_mean_accumulated', dsname='weighted', color='k', label='Weighted Avg')
 display.day_night_background('sgpmetE13.b1')
-display.axes[0].legend()
+display.axes[0, 0].legend()
 plt.show()