diff --git a/echopype/calibrate/cal_params.py b/echopype/calibrate/cal_params.py
index 040a182bef..29146de924 100644
--- a/echopype/calibrate/cal_params.py
+++ b/echopype/calibrate/cal_params.py
@@ -231,13 +231,7 @@ def _get_interp_da(
                 BB_factor.sel(channel=ch_id) if isinstance(BB_factor, xr.DataArray) else BB_factor
             )
             if isinstance(alternative, xr.DataArray):
-                # drop the redundant beam dimension if exist
-                if "beam" in alternative.coords:
-                    param.append(
-                        (alternative.sel(channel=ch_id).isel(beam=0) * BB_factor_ch).data.squeeze()
-                    )
-                else:
-                    param.append((alternative.sel(channel=ch_id) * BB_factor_ch).data.squeeze())
+                param.append((alternative.sel(channel=ch_id) * BB_factor_ch).data.squeeze())
             elif isinstance(alternative, (int, float)):
                 # expand to have ping_time dimension
                 param.append(
@@ -352,8 +346,8 @@ def get_cal_params_AZFP(beam: xr.DataArray, vend: xr.DataArray, user_dict: dict)
         if v is None:
             # Params from Sonar/Beam_group1
             if p == "equivalent_beam_angle":
-                # equivalent_beam_angle has dims: channel, ping_time, beam --> only need channel
-                out_dict[p] = beam[p].isel(ping_time=0, beam=0).drop(["ping_time", "beam"])
+                # equivalent_beam_angle has dims: channel, ping_time --> only need channel
+                out_dict[p] = beam[p].isel(ping_time=0).drop("ping_time")
 
             # Params from Vendor_specific group
             elif p in ["EL", "DS", "TVR", "VTX", "Sv_offset"]:
@@ -451,12 +445,8 @@ def _get_fs():
                 # CW: params do not require interpolation, except for impedance_transducer
                 if waveform_mode == "CW":
                     if p in PARAM_BEAM_NAME_MAP.keys():
-                        p_beam = PARAM_BEAM_NAME_MAP[p]
                         # pull from data file, these should always exist
-                        if "beam" in beam[p_beam].coords:
-                            out_dict[p] = beam[p_beam].isel(beam=0).drop("beam")
-                        else:
-                            out_dict[p] = beam[p_beam]
+                        out_dict[p] = beam[PARAM_BEAM_NAME_MAP[p]]
                     elif p == "gain_correction":
                         # pull from data file narrowband table
                         out_dict[p] = get_vend_cal_params_power(beam=beam, vend=vend, param=p)
@@ -497,7 +487,7 @@ def _get_fs():
                         )
                     elif p == "equivalent_beam_angle":
                         # scaled according to frequency ratio
-                        out_dict[p] = beam[p].isel(beam=0).drop("beam") + 20 * np.log10(
+                        out_dict[p] = beam[p] + 20 * np.log10(
                             beam["frequency_nominal"] / freq_center
                         )
                     elif p == "gain_correction":
diff --git a/echopype/calibrate/calibrate_azfp.py b/echopype/calibrate/calibrate_azfp.py
index 479b664dc9..bc2f848f03 100644
--- a/echopype/calibrate/calibrate_azfp.py
+++ b/echopype/calibrate/calibrate_azfp.py
@@ -70,9 +70,8 @@ def _cal_power_samples(self, cal_type, **kwargs):
         EL = (
             self.cal_params["EL"]
             - 2.5 / a
-            + self.echodata["Sonar/Beam_group1"]["backscatter_r"].squeeze("beam", drop=True)
-            / (26214 * a)
-        )  # eq.(5)  # has beam dim due to backscatter_r
+            + self.echodata["Sonar/Beam_group1"]["backscatter_r"] / (26214 * a)
+        )  # eq.(5)
 
         if cal_type == "Sv":
             # eq.(9)
diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py
index 0f59f2e2c0..a28b643537 100644
--- a/echopype/calibrate/calibrate_ek.py
+++ b/echopype/calibrate/calibrate_ek.py
@@ -78,7 +78,7 @@ def _cal_power_samples(self, cal_type: str) -> xr.Dataset:
             CSv = (
                 10 * np.log10(beam["transmit_power"])
                 + 2 * self.cal_params["gain_correction"]
-                + self.cal_params["equivalent_beam_angle"]  # has beam dim
+                + self.cal_params["equivalent_beam_angle"]
                 + 10
                 * np.log10(
                     wavelength**2
@@ -93,7 +93,7 @@ def _cal_power_samples(self, cal_type: str) -> xr.Dataset:
                 beam["backscatter_r"]  # has beam dim
                 + spreading_loss
                 + absorption_loss
-                - CSv  # has beam dim
+                - CSv
                 - 2 * self.cal_params["sa_correction"]
             )
             out.name = "Sv"
@@ -124,9 +124,7 @@ def _cal_power_samples(self, cal_type: str) -> xr.Dataset:
         if "time1" in out.coords:
             out = out.drop("time1")
 
-        # Squeeze out the beam dim
-        # doing it here because both out and self.cal_params["equivalent_beam_angle"] has beam dim
-        return out.squeeze("beam", drop=True)
+        return out
 
 
 class CalibrateEK60(CalibrateEK):
@@ -250,7 +248,7 @@ def __init__(
             # use true center frequency to interpolate for various cal params
             self.freq_center = (beam["frequency_start"] + beam["frequency_end"]).sel(
                 channel=self.chan_sel
-            ).isel(beam=0).drop("beam") / 2
+            ) / 2
         else:
             # use nominal channel frequency for CW pulse
             self.freq_center = beam["frequency_nominal"].sel(channel=self.chan_sel)
@@ -314,7 +312,7 @@ def _get_chan_dict(beam: xr.Dataset) -> Dict:
         if "frequency_start" in beam and "frequency_end" in beam:
             # At least some channels are BB
             # frequency_start and frequency_end are NaN for CW channels
-            freq_center = (beam["frequency_start"] + beam["frequency_end"]) / 2  # has beam dim
+            freq_center = (beam["frequency_start"] + beam["frequency_end"]) / 2
 
             return {
                 # For BB: drop channels containing CW samples (nan in freq start/end)
@@ -572,16 +570,7 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset:
         # Add env and cal parameters
         out = self._add_params_to_output(out)
 
-        # Squeeze out the beam dim
-        # out has beam dim, which came from absorption and absorption_loss
-        # self.cal_params["equivalent_beam_angle"] also has beam dim
-
-        # TODO: out should not have beam dimension at this stage
-        # once that dimension is removed from equivalent_beam_angle
-        if "beam" in out.coords:
-            return out.isel(beam=0).drop("beam")
-        else:
-            return out
+        return out
 
     def _compute_cal(self, cal_type) -> xr.Dataset:
         """
diff --git a/echopype/calibrate/range.py b/echopype/calibrate/range.py
index 8ac851f319..aa93a60616 100644
--- a/echopype/calibrate/range.py
+++ b/echopype/calibrate/range.py
@@ -183,7 +183,7 @@ def compute_range_EK(
 
     # set entries with NaN backscatter data to NaN
     if "beam" in beam["backscatter_r"].dims:
-        # Drop beam because echo_range does not have beam dimension
+        # Drop beam because echo_range should not have a beam dimension
         valid_idx = ~beam["backscatter_r"].isel(beam=0).drop("beam").isnull()
     else:
         valid_idx = ~beam["backscatter_r"].isnull()
diff --git a/echopype/consolidate/split_beam_angle.py b/echopype/consolidate/split_beam_angle.py
index 5b37b1b853..e349f399e3 100644
--- a/echopype/consolidate/split_beam_angle.py
+++ b/echopype/consolidate/split_beam_angle.py
@@ -146,12 +146,6 @@ def _e2f(angle_type: str) -> xr.Dataset:
             "from split-beam transducers!"
         )
 
-    # drop the beam dimension in theta and phi, if it exists
-    # this cannot be removed because beam dimension exists in all power data (beam dim: length=1)
-    if "beam" in theta.dims:
-        theta = theta.drop("beam").squeeze(dim="beam")
-        phi = phi.drop("beam").squeeze(dim="beam")
-
     return theta, phi
 
 
diff --git a/echopype/convert/set_groups_azfp.py b/echopype/convert/set_groups_azfp.py
index 9ae773881b..7fa404da30 100644
--- a/echopype/convert/set_groups_azfp.py
+++ b/echopype/convert/set_groups_azfp.py
@@ -19,13 +19,18 @@ class SetGroupsAZFP(SetGroupsBase):
     # these sets are applied to all Sonar/Beam_groupX groups.
 
     # Variables that need only the beam dimension added to them.
-    beam_only_names = {"backscatter_r"}
+    beam_only_names = set()
 
     # Variables that need only the ping_time dimension added to them.
-    ping_time_only_names = {"sample_interval", "transmit_duration_nominal"}
+    ping_time_only_names = {
+        "sample_interval",
+        "transmit_duration_nominal",
+        "equivalent_beam_angle",
+        "gain_correction",
+    }
 
     # Variables that need beam and ping_time dimensions added to them.
-    beam_ping_time_names = {"equivalent_beam_angle", "gain_correction"}
+    beam_ping_time_names = set()
 
     beamgroups_possible = [
         {
diff --git a/echopype/convert/set_groups_ek60.py b/echopype/convert/set_groups_ek60.py
index 7d3dcf10e1..bd99babff8 100644
--- a/echopype/convert/set_groups_ek60.py
+++ b/echopype/convert/set_groups_ek60.py
@@ -24,14 +24,16 @@ class SetGroupsEK60(SetGroupsBase):
     # in converting from v0.5.x to v0.6.0. The values within
     # these sets are applied to all Sonar/Beam_groupX groups.
 
+    # 2023-05-25 note: We have decided to remove the beam dimension from EK60,
+    # where it was added as a length-1 dimension only to more closely match
+    # the SONAR-netCDF4 v1 convention. For the time being, we are retaining the
+    # infrastructure that adds this dimension, but updating the variables lists.
+
     # Variables that need only the beam dimension added to them.
-    beam_only_names = {"backscatter_r", "angle_athwartship", "angle_alongship"}
+    beam_only_names = set()
 
     # Variables that need only the ping_time dimension added to them.
-    ping_time_only_names = {"beam_type"}
-
-    # Variables that need beam and ping_time dimensions added to them.
-    beam_ping_time_names = {
+    ping_time_only_names = {
         "beam_direction_x",
         "beam_direction_y",
         "beam_direction_z",
@@ -45,6 +47,9 @@ class SetGroupsEK60(SetGroupsBase):
         "gain_correction",
     }
 
+    # Variables that need beam and ping_time dimensions added to them.
+    beam_ping_time_names = set()
+
     beamgroups_possible = [
         {
             "name": "Beam_group1",
diff --git a/echopype/convert/set_groups_ek80.py b/echopype/convert/set_groups_ek80.py
index 308de1f1c9..de1e8f6a33 100644
--- a/echopype/convert/set_groups_ek80.py
+++ b/echopype/convert/set_groups_ek80.py
@@ -27,20 +27,11 @@ class SetGroupsEK80(SetGroupsBase):
     # these sets are applied to all Sonar/Beam_groupX groups.
 
     # Variables that need only the beam dimension added to them.
-    beam_only_names = {
-        "backscatter_r",
-        "backscatter_i",
-        "angle_athwartship",
-        "angle_alongship",
-        "frequency_start",
-        "frequency_end",
-    }
+    beam_only_names = set()
 
     # Variables that need only the ping_time dimension added to them.
-    ping_time_only_names = {"beam_type"}
-
-    # Variables that need beam and ping_time dimensions added to them.
-    beam_ping_time_names = {
+    ping_time_only_names = {
+        "beam_type",
         "beam_direction_x",
         "beam_direction_y",
         "beam_direction_z",
@@ -53,6 +44,9 @@ class SetGroupsEK80(SetGroupsBase):
         "beamwidth_twoway_athwartship",
     }
 
+    # Variables that need beam and ping_time dimensions added to them.
+    beam_ping_time_names = set()
+
     beamgroups_possible = [
         {
             "name": "Beam_group1",
diff --git a/echopype/tests/calibrate/test_cal_params.py b/echopype/tests/calibrate/test_cal_params.py
index 08e4eb05bf..aabfe2039a 100644
--- a/echopype/tests/calibrate/test_cal_params.py
+++ b/echopype/tests/calibrate/test_cal_params.py
@@ -39,11 +39,11 @@ def beam_AZFP():
     """
     beam = xr.Dataset()
     beam["equivalent_beam_angle"] = xr.DataArray(
-        [[[10, 20]]],
-        dims=["ping_time", "beam", "channel"],
-        coords={"channel": ["chA", "chB"], "ping_time": [1], "beam": [1]},
+        [[10, 20]],
+        dims=["ping_time", "channel"],
+        coords={"channel": ["chA", "chB"], "ping_time": [1]},
     )
-    return beam.transpose("channel", "ping_time", "beam")
+    return beam.transpose("channel", "ping_time")
 
 
 @pytest.fixture
@@ -84,12 +84,12 @@ def beam_EK():
         "beamwidth_twoway_alongship", "beamwidth_twoway_athwartship"
     ]:
         beam[p_name] = xr.DataArray(
-            np.array([[[123, 123, 123, 123], [456, 456, 456, 456]]]),
-            dims=["ping_time", "channel", "beam"],
-            coords={"channel": ["chA", "chB"], "ping_time": [1], "beam": [1, 2, 3, 4]},
+            np.array([[123], [456]]),
+            dims=["channel", "ping_time"],
+            coords={"channel": ["chA", "chB"], "ping_time": [1]},
         )
     beam["frequency_nominal"] = xr.DataArray([25, 55], dims=["channel"], coords={"channel": ["chA", "chB"]})
-    return beam.transpose("channel", "ping_time", "beam")
+    return beam.transpose("channel", "ping_time")
 
 
 @pytest.mark.parametrize(
@@ -262,25 +262,6 @@ def test_sanitize_user_cal_dict(sonar_type, user_dict, channel, out_dict):
                 coords={"ping_time": [1], "channel": ["chA", "chB"]}
             ),
         ),
-        #       - xr.DataArray with coordinates channel, ping_time, beam
-        (
-            xr.DataArray(
-                np.array([[np.nan, np.nan, np.nan, 4, 5, 6]]),
-                dims=["cal_channel_id", "cal_frequency"],
-                coords={"cal_channel_id": ["chB"],
-                        "cal_frequency": [10, 20, 30, 40, 50, 60]},
-            ),
-            xr.DataArray(
-                np.array([[[100, 200]]] * 4),
-                dims=["beam", "ping_time", "channel"],
-                coords={"beam": [0, 1, 2, 3], "ping_time": [1], "channel": ["chA", "chB"]},
-            ),
-            xr.DataArray(
-                [[100], [5.5]],
-                dims=["channel", "ping_time"],
-                coords={"ping_time": [1], "channel": ["chA", "chB"]}
-            ),
-        ),
         #       - xr.DataArray with coordinates channel, ping_time
         (
             xr.DataArray(
@@ -313,7 +294,6 @@ def test_sanitize_user_cal_dict(sonar_type, user_dict, channel, out_dict):
         "in_None_alt_da",
         "in_da_all_channel_out_interp",
         "in_da_some_channel_alt_scalar",
-        "in_da_some_channel_alt_da3coords",  # channel, ping_time, beam
         "in_da_some_channel_alt_da2coords",  # channel, ping_time
     ]
 )
diff --git a/echopype/tests/calibrate/test_cal_params_integration.py b/echopype/tests/calibrate/test_cal_params_integration.py
index 73bbbb303b..7a8913367a 100644
--- a/echopype/tests/calibrate/test_cal_params_integration.py
+++ b/echopype/tests/calibrate/test_cal_params_integration.py
@@ -110,9 +110,7 @@ def test_cal_params_intake_EK80_BB_complex(ek80_cal_path):
     beam = ed["Sonar/Beam_group1"].sel(channel=chan_sel)
     vend = ed["Vendor_specific"].sel(channel=chan_sel)
     freq_center = (
-        (beam["frequency_start"] + beam["frequency_end"])
-        .sel(channel=chan_sel).isel(beam=0).drop("beam") / 2
-    )
+        (beam["frequency_start"] + beam["frequency_end"]).sel(channel=chan_sel) / 2)
     cal_params_manual = ep.calibrate.cal_params.get_cal_params_EK(
         "BB", freq_center, beam, vend, {"gain_correction": gain_freq_dep}
     )
diff --git a/echopype/tests/calibrate/test_ecs_integration.py b/echopype/tests/calibrate/test_ecs_integration.py
index 9154db1427..4bf736d2ce 100644
--- a/echopype/tests/calibrate/test_ecs_integration.py
+++ b/echopype/tests/calibrate/test_ecs_integration.py
@@ -154,7 +154,7 @@ def test_ecs_intake_ek80_BB_complex(ek80_path, ecs_path):
     chan_w_BB_param = "WBT 549762-15 ES70-7C_ES"
     freq_center = (
         (beam["frequency_start"] + beam["frequency_end"]) / 2
-    ).sel(channel=chan_w_BB_param).isel(beam=0).drop_vars(["beam", "channel"])
+    ).sel(channel=chan_w_BB_param).drop_vars(["channel"])
 
     for p_name in [
         "gain_correction", "angle_offset_alongship", "angle_offset_athwartship",
diff --git a/echopype/tests/convert/test_convert_azfp.py b/echopype/tests/convert/test_convert_azfp.py
index 542587544a..075fde0fb4 100644
--- a/echopype/tests/convert/test_convert_azfp.py
+++ b/echopype/tests/convert/test_convert_azfp.py
@@ -60,7 +60,7 @@ def test_convert_azfp_01a_matlab_raw(azfp_path):
         np.array(
             [ds_matlab_output['Output'][0]['N'][fidx] for fidx in range(4)]
         ),
-        echodata["Sonar/Beam_group1"].backscatter_r.isel(beam=0).drop_vars('beam').values,
+        echodata["Sonar/Beam_group1"].backscatter_r.values,
     )
 
     # Test vendor group
@@ -128,7 +128,7 @@ def test_convert_azfp_01a_raw_echoview(azfp_path):
     echodata = open_raw(
         raw_file=azfp_01a_path, sonar_model='AZFP', xml_path=azfp_xml_path
     )
-    assert np.array_equal(test_power, echodata["Sonar/Beam_group1"].backscatter_r.isel(beam=0).drop_vars('beam')) # noqa
+    assert np.array_equal(test_power, echodata["Sonar/Beam_group1"].backscatter_r)
 
     # check convention-required variables in the Platform group
     check_platform_required_scalar_vars(echodata)
@@ -145,10 +145,10 @@ def test_convert_azfp_01a_different_ranges(azfp_path):
     )
     assert echodata["Sonar/Beam_group1"].backscatter_r.sel(channel='55030-125-1').dropna(
         'range_sample'
-    ).shape == (360, 438, 1)
+    ).shape == (360, 438)
     assert echodata["Sonar/Beam_group1"].backscatter_r.sel(channel='55030-769-4').dropna(
         'range_sample'
-    ).shape == (360, 135, 1)
+    ).shape == (360, 135)
 
     # check convention-required variables in the Platform group
     check_platform_required_scalar_vars(echodata)
diff --git a/echopype/tests/convert/test_convert_ek60.py b/echopype/tests/convert/test_convert_ek60.py
index 39a8f30c24..920775311d 100644
--- a/echopype/tests/convert/test_convert_ek60.py
+++ b/echopype/tests/convert/test_convert_ek60.py
@@ -69,7 +69,7 @@ def test_convert_ek60_matlab_raw(ek60_path):
             ds_matlab['rawData'][0]['pings'][0]['power'][0][fidx]
             for fidx in range(5)
         ],
-        echodata["Sonar/Beam_group1"].backscatter_r.isel(beam=0).transpose(
+        echodata["Sonar/Beam_group1"].backscatter_r.transpose(
             'channel', 'range_sample', 'ping_time'
         ),
         rtol=0,
@@ -82,7 +82,7 @@ def test_convert_ek60_matlab_raw(ek60_path):
                 ds_matlab['rawData'][0]['pings'][0][angle][0][fidx]
                 for fidx in range(5)
             ],
-            echodata["Sonar/Beam_group1"]['angle_' + angle].isel(beam=0).transpose(
+            echodata["Sonar/Beam_group1"]['angle_' + angle].transpose(
                 'channel', 'range_sample', 'ping_time'
             ),
         )
@@ -121,8 +121,7 @@ def test_convert_ek60_echoview_raw(ek60_path):
             echodata["Sonar/Beam_group1"].backscatter_r.isel(
                 channel=sorted_freq_ind[fidx],
                 ping_time=slice(None, 10),
-                range_sample=slice(1, None),
-                beam=0
+                range_sample=slice(1, None)
             ),
             atol=9e-6,
             rtol=atol,
diff --git a/echopype/tests/convert/test_convert_ek80.py b/echopype/tests/convert/test_convert_ek80.py
index 5354f480c6..f5ee5f6fe6 100644
--- a/echopype/tests/convert/test_convert_ek80.py
+++ b/echopype/tests/convert/test_convert_ek80.py
@@ -189,8 +189,7 @@ def test_convert_ek80_cw_power_angle_echoview(ek80_path):
         test_power = pd.read_csv(file, delimiter=';').iloc[:, 13:].values
         assert np.allclose(
             test_power,
-            echodata["Sonar/Beam_group1"].backscatter_r.sel(channel=chan,
-                                            beam='1').dropna('range_sample'),
+            echodata["Sonar/Beam_group1"].backscatter_r.sel(channel=chan).dropna('range_sample'),
             rtol=0,
             atol=1.1e-5,
         )
@@ -216,7 +215,7 @@ def test_convert_ek80_cw_power_angle_echoview(ek80_path):
         for ping_idx in df_angle['Ping_index'].value_counts().index:
             assert np.allclose(
                 df_angle.loc[df_angle['Ping_index'] == ping_idx, ' Major'],
-                major.sel(channel=chan, beam='1')
+                major.sel(channel=chan)
                 .isel(ping_time=ping_idx)
                 .dropna('range_sample'),
                 rtol=0,
@@ -224,7 +223,7 @@ def test_convert_ek80_cw_power_angle_echoview(ek80_path):
             )
             assert np.allclose(
                 df_angle.loc[df_angle['Ping_index'] == ping_idx, ' Minor'],
-                minor.sel(channel=chan, beam='1')
+                minor.sel(channel=chan)
                 .isel(ping_time=ping_idx)
                 .dropna('range_sample'),
                 rtol=0,
diff --git a/echopype/tests/echodata/test_echodata.py b/echopype/tests/echodata/test_echodata.py
index 65a16ec0c9..e78c007eb4 100644
--- a/echopype/tests/echodata/test_echodata.py
+++ b/echopype/tests/echodata/test_echodata.py
@@ -459,12 +459,14 @@ def test_nan_range_entries(range_check_files):
             echodata, env_params=None, cal_params=None, ecs_file=None, waveform_mode="BB", encode_mode="complex",
         )
         range_output = cal_obj.range_meter
+        # broadband complex data EK80 file: always need to drop "beam" dimension
         nan_locs_backscatter_r = ~echodata["Sonar/Beam_group1"].backscatter_r.isel(beam=0).drop("beam").isnull()
     else:
+        # EK60 file does not need dropping "beam" dimension
         ds_Sv = echopype.calibrate.compute_Sv(echodata)
         cal_obj = CalibrateEK60(echodata, env_params={}, cal_params=None, ecs_file=None)
         range_output = cal_obj.range_meter
-        nan_locs_backscatter_r = ~echodata["Sonar/Beam_group1"].backscatter_r.isel(beam=0).drop("beam").isnull()
+        nan_locs_backscatter_r = ~echodata["Sonar/Beam_group1"].backscatter_r.isnull()
 
     nan_locs_Sv_range = ~ds_Sv.echo_range.isnull()
     nan_locs_range = ~range_output.isnull()