diff --git a/package/CHANGELOG b/package/CHANGELOG
index 7ac8f41b18..c3002544db 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -37,6 +37,8 @@ Fixes
    DSSP by porting upstream PyDSSP 0.9.1 fix (Issue #4913)
 
 Enhancements
+ * Enables parallelization for analysis.hydrogenbonds.wbridge_analysis.WaterBridgeAnalysis
+   (Issue #4666, PR #5151)
  * Adds support for parsing `.tpr` files produced by GROMACS 2026.0
  * Enables parallelization for analysis.diffusionmap.DistanceMatrix
    (Issue #4679, PR #4745)
diff --git a/package/MDAnalysis/analysis/hydrogenbonds/wbridge_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/wbridge_analysis.py
index 027c0b7125..aba9bd152b 100644
--- a/package/MDAnalysis/analysis/hydrogenbonds/wbridge_analysis.py
+++ b/package/MDAnalysis/analysis/hydrogenbonds/wbridge_analysis.py
@@ -716,7 +716,7 @@ def analysis(current, output, u, **kwargs):
 from MDAnalysis.lib.distances import calc_angles, capped_distance
 from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch
 
-from ..base import AnalysisBase
+from ..base import AnalysisBase, ResultsGroup
 
 logger = logging.getLogger("MDAnalysis.analysis.WaterBridgeAnalysis")
 
@@ -731,6 +731,10 @@ class WaterBridgeAnalysis(AnalysisBase):
 
     .. versionadded:: 0.17.0
 
+    .. versionchanged:: 2.11.0
+       Enabled **parallel execution** with the ``multiprocessing`` and ``dask``
+       backends; use the new method :meth:`get_supported_backends` to see all
+       supported backends.
     """
 
     # use tuple(set()) here so that one can just copy&paste names from the
@@ -804,6 +808,16 @@ class WaterBridgeAnalysis(AnalysisBase):
         lambda: 1.5, N=1.31, O=1.31, P=1.58, S=1.55  # default value
     )  # noqa: E741
 
+    _analysis_algorithm_is_parallelizable = True
+
+    @classmethod
+    def get_supported_backends(cls):
+        return (
+            "serial",
+            "multiprocessing",
+            "dask",
+        )
+
     def __init__(
         self,
         universe,
@@ -1014,7 +1028,6 @@ def __init__(
         # final result accessed as self.results.network
         self.results.network = []
         self.results.timeseries = None
-        self.timesteps = None  # time for each frame
 
         self._log_parameters()
 
@@ -1301,7 +1314,6 @@ def _prepare(self):
 
         self._update_selection()
 
-        self.timesteps = []
         if len(self._s1) and len(self._s2):
             self._update_water_selection()
         else:
@@ -1396,7 +1408,6 @@ def _donor2acceptor(self, donors, h_donors, acceptor):
         return result
 
     def _single_frame(self):
-        self.timesteps.append(self._ts.time)
         self.box = self.u.dimensions if self.pbc else None
 
         if self.update_selection:
@@ -1953,7 +1964,8 @@ def count_by_time(self, analysis_func=None, **kwargs):
             analysis_func = self._count_by_time_analysis
         if self.results.network:
             result = []
-            for time, frame in zip(self.timesteps, self.results.network):
+
+            for time, frame in zip(self.times, self.results.network):
                 result_dict = defaultdict(int)
                 self._traverse_water_network(
                     frame,
@@ -2016,11 +2028,8 @@ def timesteps_by_type(self, analysis_func=None, **kwargs):
 
         if self.results.network:
             result = defaultdict(list)
-            if self.timesteps is None:
-                timesteps = range(len(self.results.network))
-            else:
-                timesteps = self.timesteps
-            for time, frame in zip(timesteps, self.results.network):
+
+            for time, frame in zip(self.times, self.results.network):
                 self._traverse_water_network(
                     frame,
                     [],
@@ -2120,7 +2129,8 @@ def generate_table(self, output_format=None):
         # standard array, like this:
         out = np.empty((num_records,), dtype=dtype)
         cursor = 0  # current row
-        for t, hframe in zip(self.timesteps, timeseries):
+
+        for t, hframe in zip(self.times, timeseries):
             for (
                 donor_index,
                 acceptor_index,
@@ -2152,6 +2162,14 @@ def generate_table(self, output_format=None):
     def _conclude(self):
         self.results.timeseries = self._generate_timeseries()
 
+    def _get_aggregator(self):
+        return ResultsGroup(
+            lookup={
+                "timeseries": ResultsGroup.ndarray_hstack,
+                "network": ResultsGroup.ndarray_hstack,
+            }
+        )
+
     @property
     def network(self):
         wmsg = (
diff --git a/testsuite/MDAnalysisTests/analysis/conftest.py b/testsuite/MDAnalysisTests/analysis/conftest.py
index 7ffe99ef32..067ee5702c 100644
--- a/testsuite/MDAnalysisTests/analysis/conftest.py
+++ b/testsuite/MDAnalysisTests/analysis/conftest.py
@@ -15,6 +15,9 @@
 from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import (
     HydrogenBondAnalysis,
 )
+from MDAnalysis.analysis.hydrogenbonds.wbridge_analysis import (
+    WaterBridgeAnalysis,
+)
 from MDAnalysis.analysis.nucleicacids import NucPairDist
 from MDAnalysis.analysis.contacts import Contacts
 from MDAnalysis.analysis.density import DensityAnalysis
@@ -217,3 +220,11 @@ def client_InterRDF_s(request):
 @pytest.fixture(scope="module", params=params_for_cls(DistanceMatrix))
 def client_DistanceMatrix(request):
     return request.param
+
+
+# MDAnalysis.analysis.hydrogenbonds.wbridge_analysis
+
+
+@pytest.fixture(scope="module", params=params_for_cls(WaterBridgeAnalysis))
+def client_WaterBridgeAnalysis(request):
+    return request.param
diff --git a/testsuite/MDAnalysisTests/analysis/test_wbridge.py b/testsuite/MDAnalysisTests/analysis/test_wbridge.py
index 2b53d2f1d3..97ae26889a 100644
--- a/testsuite/MDAnalysisTests/analysis/test_wbridge.py
+++ b/testsuite/MDAnalysisTests/analysis/test_wbridge.py
@@ -1,4 +1,4 @@
-from io import StringIO
+import numpy as np
 from collections import defaultdict
 
 from numpy.testing import (
@@ -7,273 +7,43 @@
 )
 import pytest
 
+from MDAnalysis.lib.util import NamedStream
+
 import MDAnalysis
 from MDAnalysis.analysis.hydrogenbonds.wbridge_analysis import (
     WaterBridgeAnalysis,
 )
 
 
-class TestWaterBridgeAnalysis(object):
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_empty():
-        """A universe with no hydrogen bonds"""
-        grofile = """Test gro file
-5
-    1ALA      N    1   0.000   0.000   0.000
-    1ALA      H    2   0.100   0.000   0.000
-    2SOL     OW    3   3.000   0.000   0.000
-    4ALA      H    4   0.500   0.000   0.000
-    4ALA      N    5   0.600   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_DA():
-        """A universe with one hydrogen bond acceptor bonding to a hydrogen bond
-        donor"""
-        grofile = """Test gro file
-3
-    1ALA      N    1   0.000   0.000   0.000
-    1ALA      H    2   0.100   0.000   0.000
-    4ALA      O    3   0.300   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_DA_PBC():
-        """A universe with one hydrogen bond acceptor bonding to a hydrogen bond
-        donor but in a PBC condition"""
-        grofile = """Test gro file
-3
-    1ALA      N    1   0.800   0.000   0.000
-    1ALA      H    2   0.900   0.000   0.000
-    4ALA      O    3   0.100   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_AD():
-        """A universe with one hydrogen bond donor bonding to a hydrogen bond
-        acceptor"""
-        grofile = """Test gro file
-3
-    1ALA      O    1   0.000   0.000   0.000
-    4ALA      H    2   0.200   0.000   0.000
-    4ALA      N    3   0.300   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_loop():
-        """A universe with one hydrogen bond acceptor bonding to a water which
-        bonds back to the first hydrogen bond  acceptor and thus form a loop"""
-        grofile = """Test gro file
-5
-    1ALA      O    1   0.000   0.001   0.000
-    2SOL     OW    2   0.300   0.001   0.000
-    2SOL    HW1    3   0.200   0.002   0.000
-    2SOL    HW2    4   0.200   0.000   0.000
-    4ALA      O    5   0.600   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_DWA():
-        """A universe with one hydrogen bond donor bonding to a hydrogen bond
-        acceptor through a water"""
-        grofile = """Test gro file
-5
-    1ALA      N    1   0.000   0.000   0.000
-    1ALA      H    2   0.100   0.000   0.000
-    2SOL     OW    3   0.300   0.000   0.000
-    2SOL    HW2    4   0.400   0.000   0.000
-    4ALA      O    5   0.600   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_DWD():
-        """A universe with one hydrogen bond donor bonding to a hydrogen bond
-        donor through a water"""
-        grofile = """Test gro file
-5
-    1ALA      N    1   0.000   0.000   0.000
-    1ALA      H    2   0.100   0.000   0.000
-    2SOL     OW    3   0.300   0.000   0.000
-    4ALA      H    4   0.500   0.000   0.000
-    4ALA      N    5   0.600   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_AWA():
-        """A universe with two hydrogen bond acceptor are joined by a water"""
-        grofile = """Test gro file
-5
-    1ALA      O    1   0.000   0.000   0.000
-    2SOL     OW    2   0.300   0.000   0.000
-    2SOL    HW1    3   0.200   0.000   0.000
-    2SOL    HW2    4   0.400   0.000   0.000
-    4ALA      O    5   0.600   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_AWD():
-        """A universe with one hydrogen bond acceptor bonding to a hydrogen
-        bond donor through a water"""
-        grofile = """Test gro file
-5
-    1ALA      O    1   0.000   0.000   0.000
-    2SOL     OW    2   0.300   0.000   0.000
-    2SOL    HW1    3   0.200   0.000   0.000
-    4ALA      H    4   0.500   0.000   0.000
-    4ALA      N    5   0.600   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_AWWA():
-        """A universe with one hydrogen bond acceptor bonding to a hydrogen bond
-        acceptor through two waters"""
-        grofile = """Test gro file
-7
-    1ALA      O    1   0.000   0.000   0.000
-    2SOL     OW    2   0.300   0.000   0.000
-    2SOL    HW1    3   0.200   0.000   0.000
-    2SOL    HW2    4   0.400   0.000   0.000
-    3SOL     OW    5   0.600   0.000   0.000
-    3SOL    HW1    6   0.700   0.000   0.000
-    4ALA      O    7   0.900   0.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_AWWWA():
-        """A universe with one hydrogen bond acceptor bonding to a hydrogen bond
-        acceptor through three waters"""
-        grofile = """Test gro file
-9
-    1ALA      O    1   0.000   0.000   0.000
-    2SOL     OW    2   0.300   0.000   0.000
-    2SOL    HW1    3   0.200   0.000   0.000
-    2SOL    HW2    4   0.400   0.000   0.000
-    3SOL     OW    5   0.600   0.000   0.000
-    3SOL    HW1    6   0.700   0.000   0.000
-    4SOL     OW    7   0.900   0.000   0.000
-    4SOL    HW1    8   1.000   0.000   0.000
-    5ALA      O    9   1.200   0.000   0.000
-  10.0   10.0   10.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_AWWWWA():
-        """A universe with one hydrogen bond acceptor bonding to a hydrogen bond
-        acceptor through three waters"""
-        grofile = """Test gro file
-11
-    1ALA      O    1   0.000   0.000   0.000
-    2SOL     OW    2   0.300   0.000   0.000
-    2SOL    HW1    3   0.200   0.000   0.000
-    2SOL    HW2    4   0.400   0.000   0.000
-    3SOL     OW    5   0.600   0.000   0.000
-    3SOL    HW1    6   0.700   0.000   0.000
-    4SOL     OW    7   0.900   0.000   0.000
-    4SOL    HW1    8   1.000   0.000   0.000
-    5SOL     OW    9   1.200   0.000   0.000
-    5SOL    HW1   10   1.300   0.000   0.000
-    6ALA      O   11   1.400   0.000   0.000
-  10.0   10.0   10.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
-
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_branch():
-        """A universe with one hydrogen bond acceptor bonding to two hydrogen
-        bond acceptor in selection 2"""
-        grofile = """Test gro file
-9
-    1ALA      O    1   0.000   0.000   0.000
-    2SOL     OW    2   0.300   0.000   0.000
-    2SOL    HW1    3   0.200   0.000   0.000
-    2SOL    HW2    4   0.400   0.000   0.000
-    3SOL     OW    5   0.600   0.000   0.000
-    3SOL    HW1    6   0.700   0.000   0.000
-    3SOL    HW2    7   0.600   0.100   0.000
-    4ALA      O    8   0.900   0.000   0.000
-    5ALA      O    9   0.600   0.300   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
+# Note: Currently the datafiles are added as a go-around for the
+# incompatibility with StringIO with a fix in the future. Issue #5221
+from MDAnalysisTests.datafiles import (
+    WB_AD,
+    WB_AWA,
+    WB_AWA_AWWA,
+    WB_AWD,
+    WB_AWWA,
+    WB_AWWWA,
+    WB_AWWWWA,
+    WB_BRANCH,
+    WB_DA,
+    WB_DA_PBC,
+    WB_DWA,
+    WB_DWD,
+    WB_EMPTY,
+    WB_LOOP,
+    WB_DUPLICATE_WATER,
+    WB_MULTIFRAME_GRO,
+    WB_MULTIFRAME_DCD,
+)
 
-    @staticmethod
-    @pytest.fixture(scope="class")
-    def universe_AWA_AWWA():
-        """A universe with one hydrogen bond acceptors are bonded through one or
-        two water"""
-        grofile = """Test gro file
-12
-    1ALA      O    1   0.000   0.000   0.000
-    2SOL     OW    2   0.300   0.000   0.000
-    2SOL    HW1    3   0.200   0.000   0.000
-    2SOL    HW2    4   0.400   0.000   0.000
-    4ALA      O    5   0.600   0.000   0.000
-    5ALA      O    6   0.000   1.000   0.000
-    6SOL     OW    7   0.300   1.000   0.000
-    6SOL    HW1    8   0.200   1.000   0.000
-    6SOL    HW2    9   0.400   1.000   0.000
-    7SOL     OW   10   0.600   1.000   0.000
-    7SOL    HW1   11   0.700   1.000   0.000
-    8ALA      O   12   0.900   1.000   0.000
-  1.0   1.0   1.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
-        return u
 
+class TestWaterBridgeAnalysis(object):
     @staticmethod
     @pytest.fixture(scope="class")
     def wb_multiframe():
         """A water bridge object with multipley frames"""
-        grofile = """Test gro file
- 13
-    1ALA      O    1   0.000   0.000   0.000
-    1ALA      H    2   0.000   0.000   0.000
-    2SOL     OW    3   0.300   0.000   0.000
-    2SOL    HW1    4   0.200   0.000   0.000
-    2SOL    HW2    5   0.400   0.000   0.000
-    3SOL     OW    6   0.600   0.000   0.000
-    3SOL    HW1    7   0.700   0.000   0.000
-    4SOL     OW    8   0.900   0.000   0.000
-    4SOL    HW1    9   1.000   0.000   0.000
-    5SOL     OW   10   1.200   0.000   0.000
-    5SOL    HW1   11   1.300   0.000   0.000
-    6ALA      H   12   1.400   0.000   0.000
-    6ALA      O   13   1.400   0.000   0.000
-   10.0   10.0   10.0"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
+        u = MDAnalysis.Universe(WB_MULTIFRAME_GRO, WB_MULTIFRAME_DCD)
         wb = WaterBridgeAnalysis(
             u, "protein and (resid 1)", "protein and (resid 4)", order=4
         )
@@ -293,12 +63,13 @@ def wb_multiframe():
                 }
             }
         )
-        wb.timesteps = range(len(wb.results.network))
+        wb.times = np.arange(len(wb.results.network))
         return wb
 
-    def test_nodata(self, universe_DA):
+    def test_nodata(self):
         """Test if the funtions can run when there is no data.
         This is achieved by not runing the run() first."""
+        universe_DA = MDAnalysis.Universe(WB_DA)
         wb = WaterBridgeAnalysis(
             universe_DA,
             "protein and (resid 1)",
@@ -310,8 +81,9 @@ def test_nodata(self, universe_DA):
         assert_equal(wb.count_by_time(), None)
         assert_equal(wb.count_by_type(), None)
 
-    def test_selection_type_error(self, universe_DA):
+    def test_selection_type_error(self):
         """Test the case when the wrong selection1_type is given"""
+        universe_DA = MDAnalysis.Universe(WB_DA)
         try:
             wb = WaterBridgeAnalysis(
                 universe_DA,
@@ -325,8 +97,9 @@ def test_selection_type_error(self, universe_DA):
         else:
             raise pytest.fail("selection_type aaa should rasie error")
 
-    def test_distance_type_error(self, universe_DA):
+    def test_distance_type_error(self):
         """Test the case when the wrong selection1_type is given"""
+        universe_DA = MDAnalysis.Universe(WB_DA)
         with pytest.raises(
             ValueError,
             match="Only 'hydrogen' and 'heavy' are allowed for option `distance_type'",
@@ -340,8 +113,9 @@ def test_distance_type_error(self, universe_DA):
                 distance_type="aaa",
             )
 
-    def test_selection2_type_error(self, universe_DA):
+    def test_selection2_type_error(self):
         """Test the case when the wrong selection1_type is given"""
+        universe_DA = MDAnalysis.Universe(WB_DA)
         with pytest.raises(
             ValueError, match="`selection2_type` is not a keyword argument."
         ):
@@ -354,30 +128,33 @@ def test_selection2_type_error(self, universe_DA):
                 selection2_type="aaa",
             )
 
-    def test_empty_selection(self, universe_DA):
+    def test_empty_selection(self, client_WaterBridgeAnalysis):
         """Test the case when selection yields empty result"""
+        universe_DA = MDAnalysis.Universe(WB_DA)
         wb = WaterBridgeAnalysis(
             universe_DA,
             "protein and (resid 9)",
             "protein and (resid 10)",
             order=0,
         )
-        wb.run()
+        wb.run(**client_WaterBridgeAnalysis)
         assert wb.results.network == [{}]
 
-    def test_loop(self, universe_loop):
+    def test_loop(self, client_WaterBridgeAnalysis):
         """Test if loop can be handled correctly"""
+        universe_loop = MDAnalysis.Universe(WB_LOOP)
         wb = WaterBridgeAnalysis(
             universe_loop,
             "protein and (resid 1)",
             "protein and (resid 1 or resid 4)",
         )
-        wb.run()
+        wb.run(**client_WaterBridgeAnalysis)
         assert_equal(len(wb.results.network[0].keys()), 2)
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_donor_accepter(self, universe_DA, distance_type):
+    def test_donor_accepter(self, distance_type, client_WaterBridgeAnalysis):
         """Test zeroth order donor to acceptor hydrogen bonding"""
+        universe_DA = MDAnalysis.Universe(WB_DA)
         wb = WaterBridgeAnalysis(
             universe_DA,
             "protein and (resid 1)",
@@ -387,13 +164,16 @@ def test_donor_accepter(self, universe_DA, distance_type):
             debug=True,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (1, 0, 2, None))
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_donor_accepter_pbc(self, universe_DA_PBC, distance_type):
+    def test_donor_accepter_pbc(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test zeroth order donor to acceptor hydrogen bonding in PBC conditions"""
+        universe_DA_PBC = MDAnalysis.Universe(WB_DA_PBC)
         wb = WaterBridgeAnalysis(
             universe_DA_PBC,
             "protein and (resid 1)",
@@ -402,13 +182,14 @@ def test_donor_accepter_pbc(self, universe_DA_PBC, distance_type):
             pbc=True,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (1, 0, 2, None))
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_accepter_donor(self, universe_AD, distance_type):
+    def test_accepter_donor(self, distance_type, client_WaterBridgeAnalysis):
         """Test zeroth order acceptor to donor hydrogen bonding"""
+        universe_AD = MDAnalysis.Universe(WB_AD)
         wb = WaterBridgeAnalysis(
             universe_AD,
             "protein and (resid 1)",
@@ -416,21 +197,24 @@ def test_accepter_donor(self, universe_AD, distance_type):
             order=0,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 1, 2))
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_acceptor_water_accepter(self, universe_AWA, distance_type):
+    def test_acceptor_water_accepter(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test case where the hydrogen bond acceptor from selection 1 form
         water bridge with hydrogen bond acceptor from selection 2"""
+        universe_AWA = MDAnalysis.Universe(WB_AWA)
         wb = WaterBridgeAnalysis(
             universe_AWA,
             "protein and (resid 1)",
             "protein and (resid 4)",
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -438,16 +222,19 @@ def test_acceptor_water_accepter(self, universe_AWA, distance_type):
         assert_equal(second[list(second.keys())[0]], None)
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_donor_water_accepter(self, universe_DWA, distance_type):
+    def test_donor_water_accepter(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test case where the hydrogen bond donor from selection 1 form
         water bridge with hydrogen bond acceptor from selection 2"""
+        universe_DWA = MDAnalysis.Universe(WB_DWA)
         wb = WaterBridgeAnalysis(
             universe_DWA,
             "protein and (resid 1)",
             "protein and (resid 4)",
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (1, 0, 2, None))
         second = network[list(network.keys())[0]]
@@ -455,16 +242,19 @@ def test_donor_water_accepter(self, universe_DWA, distance_type):
         assert_equal(second[list(second.keys())[0]], None)
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_acceptor_water_donor(self, universe_AWD, distance_type):
+    def test_acceptor_water_donor(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test case where the hydrogen bond acceptor from selection 1 form
         water bridge with hydrogen bond donor from selection 2"""
+        universe_AWD = MDAnalysis.Universe(WB_AWD)
         wb = WaterBridgeAnalysis(
             universe_AWD,
             "protein and (resid 1)",
             "protein and (resid 4)",
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -472,52 +262,60 @@ def test_acceptor_water_donor(self, universe_AWD, distance_type):
         assert_equal(second[list(second.keys())[0]], None)
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_donor_water_donor(self, universe_DWD, distance_type):
+    def test_donor_water_donor(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test case where the hydrogen bond donor from selection 1 form
         water bridge with hydrogen bond donor from selection 2"""
+        universe_DWD = MDAnalysis.Universe(WB_DWD)
         wb = WaterBridgeAnalysis(
             universe_DWD,
             "protein and (resid 1)",
             "protein and (resid 4)",
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (1, 0, 2, None))
         second = network[list(network.keys())[0]]
         assert_equal(list(second.keys())[0][:4], (2, None, 3, 4))
         assert_equal(second[list(second.keys())[0]], None)
 
-    def test_empty(self, universe_empty):
+    def test_empty(self, client_WaterBridgeAnalysis):
         """Test case where no water bridge exists"""
+        universe_empty = MDAnalysis.Universe(WB_EMPTY)
         wb = WaterBridgeAnalysis(universe_empty, "protein", "protein")
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         assert_equal(wb.results.network[0], defaultdict(dict))
 
-    def test_same_selection(self, universe_DWA):
+    def test_same_selection(self, client_WaterBridgeAnalysis):
         """
         This test tests that if the selection 1 and selection 2 are both protein.
         However, the protein only forms one hydrogen bond with the water.
         This entry won't be included.
         """
+        universe_DWA = MDAnalysis.Universe(WB_DWA)
         wb = WaterBridgeAnalysis(
             universe_DWA, "protein and resid 1", "protein and resid 1"
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         assert_equal(wb.results.network[0], defaultdict(dict))
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_acceptor_2water_accepter(self, universe_AWWA, distance_type):
+    def test_acceptor_2water_accepter(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test case where the hydrogen bond acceptor from selection 1 form second order
         water bridge with hydrogen bond acceptor from selection 2"""
         # test first order
+        universe_AWWA = MDAnalysis.Universe(WB_AWWA)
         wb = WaterBridgeAnalysis(
             universe_AWWA,
             "protein and (resid 1)",
             "protein and (resid 4)",
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         assert_equal(wb.results.network[0], defaultdict(dict))
         # test second order
         wb = WaterBridgeAnalysis(
@@ -527,7 +325,7 @@ def test_acceptor_2water_accepter(self, universe_AWWA, distance_type):
             order=2,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -543,7 +341,7 @@ def test_acceptor_2water_accepter(self, universe_AWWA, distance_type):
             order=3,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -553,9 +351,12 @@ def test_acceptor_2water_accepter(self, universe_AWWA, distance_type):
         assert_equal(third[list(third.keys())[0]], None)
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_acceptor_3water_accepter(self, universe_AWWWA, distance_type):
+    def test_acceptor_3water_accepter(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test case where the hydrogen bond acceptor from selection 1 form third order
         water bridge with hydrogen bond acceptor from selection 2"""
+        universe_AWWWA = MDAnalysis.Universe(WB_AWWWA)
         wb = WaterBridgeAnalysis(
             universe_AWWWA,
             "protein and (resid 1)",
@@ -563,7 +364,7 @@ def test_acceptor_3water_accepter(self, universe_AWWWA, distance_type):
             order=2,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         assert_equal(wb.results.network[0], defaultdict(dict))
 
         wb = WaterBridgeAnalysis(
@@ -573,7 +374,7 @@ def test_acceptor_3water_accepter(self, universe_AWWWA, distance_type):
             order=3,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -591,7 +392,7 @@ def test_acceptor_3water_accepter(self, universe_AWWWA, distance_type):
             order=4,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -603,9 +404,12 @@ def test_acceptor_3water_accepter(self, universe_AWWWA, distance_type):
         assert_equal(fourth[list(fourth.keys())[0]], None)
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_acceptor_4water_accepter(self, universe_AWWWWA, distance_type):
+    def test_acceptor_4water_accepter(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test case where the hydrogen bond acceptor from selection 1 form fourth order
         water bridge with hydrogen bond acceptor from selection 2"""
+        universe_AWWWWA = MDAnalysis.Universe(WB_AWWWWA)
         wb = WaterBridgeAnalysis(
             universe_AWWWWA,
             "protein and (resid 1)",
@@ -613,7 +417,7 @@ def test_acceptor_4water_accepter(self, universe_AWWWWA, distance_type):
             order=3,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         assert_equal(wb.results.network[0], defaultdict(dict))
 
         wb = WaterBridgeAnalysis(
@@ -623,7 +427,7 @@ def test_acceptor_4water_accepter(self, universe_AWWWWA, distance_type):
             order=4,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -643,7 +447,7 @@ def test_acceptor_4water_accepter(self, universe_AWWWWA, distance_type):
             order=5,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -657,10 +461,13 @@ def test_acceptor_4water_accepter(self, universe_AWWWWA, distance_type):
         assert_equal(fifth[list(fifth.keys())[0]], None)
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_acceptor_22water_accepter(self, universe_branch, distance_type):
+    def test_acceptor_22water_accepter(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test case where the hydrogen bond acceptor from selection 1 form a second order
         water bridge with hydrogen bond acceptor from selection 2
         and the last water is linked to two residues in selection 2"""
+        universe_branch = MDAnalysis.Universe(WB_BRANCH)
         wb = WaterBridgeAnalysis(
             universe_branch,
             "protein and (resid 1)",
@@ -668,7 +475,7 @@ def test_acceptor_22water_accepter(self, universe_branch, distance_type):
             order=2,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -679,8 +486,9 @@ def test_acceptor_22water_accepter(self, universe_branch, distance_type):
             sorted([key[:4] for key in list(third.keys())]),
         )
 
-    def test_timeseries_wba(self, universe_branch):
+    def test_timeseries_wba(self, client_WaterBridgeAnalysis):
         """Test if the time series data is correctly generated in water bridge analysis format"""
+        universe_branch = MDAnalysis.Universe(WB_BRANCH)
         wb = WaterBridgeAnalysis(
             universe_branch,
             "protein and (resid 1)",
@@ -688,7 +496,7 @@ def test_timeseries_wba(self, universe_branch):
             order=2,
         )
         wb.output_format = "sele1_sele2"
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         timeseries = sorted(wb.results.timeseries[0])
 
         assert_equal(
@@ -704,8 +512,9 @@ def test_timeseries_wba(self, universe_branch):
             timeseries[3][:4], (6, 8, ("SOL", 3, "HW2"), ("ALA", 5, "O"))
         )
 
-    def test_timeseries_hba(self, universe_branch):
+    def test_timeseries_hba(self, client_WaterBridgeAnalysis):
         """Test if the time series data is correctly generated in hydrogen bond analysis format"""
+        universe_branch = MDAnalysis.Universe(WB_BRANCH)
         wb = WaterBridgeAnalysis(
             universe_branch,
             "protein and (resid 1)",
@@ -713,7 +522,7 @@ def test_timeseries_hba(self, universe_branch):
             order=2,
         )
         wb.output_format = "donor_acceptor"
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         timeseries = sorted(wb.results.timeseries[0])
 
         assert_equal(
@@ -730,8 +539,11 @@ def test_timeseries_hba(self, universe_branch):
         )
 
     @pytest.mark.parametrize("distance_type", ["hydrogen", "heavy"])
-    def test_acceptor_12water_accepter(self, universe_AWA_AWWA, distance_type):
+    def test_acceptor_12water_accepter(
+        self, distance_type, client_WaterBridgeAnalysis
+    ):
         """Test of independent first order and second can be recognised correctely"""
+        universe_AWA_AWWA = MDAnalysis.Universe(WB_AWA_AWWA)
         wb = WaterBridgeAnalysis(
             universe_AWA_AWWA,
             "protein and (resid 1 or resid 5)",
@@ -739,7 +551,7 @@ def test_acceptor_12water_accepter(self, universe_AWA_AWWA, distance_type):
             order=1,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(list(network.keys())[0][:4], (0, None, 2, 1))
         second = network[list(network.keys())[0]]
@@ -753,36 +565,38 @@ def test_acceptor_12water_accepter(self, universe_AWA_AWWA, distance_type):
             order=2,
             distance_type=distance_type,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         network = wb.results.network[0]
         assert_equal(
             [(0, None, 2, 1), (5, None, 7, 6)],
             sorted([key[:4] for key in list(network.keys())]),
         )
 
-    def test_count_by_type_single_link(self, universe_DWA):
+    def test_count_by_type_single_link(self, client_WaterBridgeAnalysis):
         """
         This test tests the simplest water bridge to see if count_by_type() works.
         """
+        universe_DWA = MDAnalysis.Universe(WB_DWA)
         wb = WaterBridgeAnalysis(
             universe_DWA, "protein and (resid 1)", "protein and (resid 4)"
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         assert_equal(
             wb.count_by_type(), [(1, 4, "ALA", 1, "H", "ALA", 4, "O", 1.0)]
         )
 
-    def test_count_by_type_multiple_link(self, universe_AWA_AWWA):
+    def test_count_by_type_multiple_link(self, client_WaterBridgeAnalysis):
         """
         This test tests if count_by_type() can give the correct result for more than 1 links.
         """
+        universe_AWA_AWWA = MDAnalysis.Universe(WB_AWA_AWWA)
         wb = WaterBridgeAnalysis(
             universe_AWA_AWWA,
             "protein and (resid 1 or resid 5)",
             "protein and (resid 4 or resid 8)",
             order=2,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
         assert_equal(
             sorted(wb.count_by_type()),
             [
@@ -936,19 +750,20 @@ def test_count_by_time(self, wb_multiframe):
             wb_multiframe.count_by_time(), [(0, 1), (1, 1), (2, 1), (3, 1)]
         )
 
-    def test_count_by_time_weight(self, universe_AWA_AWWA):
+    def test_count_by_time_weight(self, client_WaterBridgeAnalysis):
         """
         This test tests if modyfing the analysis_func allows the weight to be changed
         in count_by_type().
         :return:
         """
+        universe_AWA_AWWA = MDAnalysis.Universe(WB_AWA_AWWA)
         wb = WaterBridgeAnalysis(
             universe_AWA_AWWA,
             "protein and (resid 1 or resid 5)",
             "protein and (resid 4 or resid 8)",
             order=2,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
 
         def analysis(current, output, u):
             sele1_index, sele1_heavy_index, atom2, heavy_atom2, dist, angle = (
@@ -979,18 +794,19 @@ def analysis(current, output, u):
             ],
         )
 
-    def test_count_by_time_empty(self, universe_AWA_AWWA):
+    def test_count_by_time_empty(self, client_WaterBridgeAnalysis):
         """
         See if count_by_time() can handle zero well.
         :return:
         """
+        universe_AWA_AWWA = MDAnalysis.Universe(WB_AWA_AWWA)
         wb = WaterBridgeAnalysis(
             universe_AWA_AWWA,
             "protein and (resid 1 or resid 5)",
             "protein and (resid 4 or resid 8)",
             order=2,
         )
-        wb.run(verbose=False)
+        wb.run(**client_WaterBridgeAnalysis, verbose=False)
 
         def analysis(current, output, u):
             pass
@@ -1026,38 +842,23 @@ def test_timesteps_by_type(self, wb_multiframe):
             timesteps[3], [1, 12, "ALA", 1, "H", "ALA", 6, "O", 0, 2]
         )
 
-    def test_duplicate_water(self):
-        """A case #3119 where
-        Acceptor···H−O···H-Donor
-                     |
-                     H···O-H
-        will be recognised as 3rd order water bridge.
-        """
-        grofile = """Test gro file
-    7
-    1LEU      O    1   1.876   0.810   1.354
-  117SOL    HW1    2   1.853   0.831   1.162
-  117SOL     OW    3   1.877   0.890   1.081
-  117SOL    HW2    4   1.908   0.828   1.007
-  135SOL     OW    5   1.924   0.713   0.845
-    1LEU      H    6   1.997   0.991   1.194
-    1LEU      N    7   2.041   1.030   1.274
-   2.22092   2.22092   2.22092"""
-        u = MDAnalysis.Universe(StringIO(grofile), format="gro")
+    def test_duplicate_water(self, client_WaterBridgeAnalysis):
+        u = MDAnalysis.Universe(WB_DUPLICATE_WATER)
         wb = WaterBridgeAnalysis(
             u, "resname LEU and name O", "resname LEU and name N H", order=4
         )
-        wb.run()
+        wb.run(**client_WaterBridgeAnalysis)
         assert len(wb.results.timeseries[0]) == 2
 
-    def test_warn_results_deprecated(self, universe_DA):
+    def test_warn_results_deprecated(self, client_WaterBridgeAnalysis):
+        universe_DA = MDAnalysis.Universe(WB_DA)
         wb = WaterBridgeAnalysis(
             universe_DA,
             "protein and (resid 9)",
             "protein and (resid 10)",
             order=0,
         )
-        wb.run()
+        wb.run(**client_WaterBridgeAnalysis)
 
         wmsg = "The `network` attribute was deprecated in MDAnalysis 2.0.0"
         with pytest.warns(DeprecationWarning, match=wmsg):
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_ad.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_ad.gro
new file mode 100644
index 0000000000..62058bb253
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_ad.gro
@@ -0,0 +1,6 @@
+Test gro file
+3
+    1ALA      O    1   0.000   0.000   0.000
+    4ALA      H    2   0.200   0.000   0.000
+    4ALA      N    3   0.300   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_awa.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_awa.gro
new file mode 100644
index 0000000000..cc1cd188bc
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_awa.gro
@@ -0,0 +1,8 @@
+Test gro file
+5
+    1ALA      O    1   0.000   0.000   0.000
+    2SOL     OW    2   0.300   0.000   0.000
+    2SOL    HW1    3   0.200   0.000   0.000
+    2SOL    HW2    4   0.400   0.000   0.000
+    4ALA      O    5   0.600   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_awa_awwa.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_awa_awwa.gro
new file mode 100644
index 0000000000..d40311bd07
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_awa_awwa.gro
@@ -0,0 +1,15 @@
+Test gro file
+12
+    1ALA      O    1   0.000   0.000   0.000
+    2SOL     OW    2   0.300   0.000   0.000
+    2SOL    HW1    3   0.200   0.000   0.000
+    2SOL    HW2    4   0.400   0.000   0.000
+    4ALA      O    5   0.600   0.000   0.000
+    5ALA      O    6   0.000   1.000   0.000
+    6SOL     OW    7   0.300   1.000   0.000
+    6SOL    HW1    8   0.200   1.000   0.000
+    6SOL    HW2    9   0.400   1.000   0.000
+    7SOL     OW   10   0.600   1.000   0.000
+    7SOL    HW1   11   0.700   1.000   0.000
+    8ALA      O   12   0.900   1.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_awd.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_awd.gro
new file mode 100644
index 0000000000..4e2488d622
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_awd.gro
@@ -0,0 +1,8 @@
+Test gro file
+5
+    1ALA      O    1   0.000   0.000   0.000
+    2SOL     OW    2   0.300   0.000   0.000
+    2SOL    HW1    3   0.200   0.000   0.000
+    4ALA      H    4   0.500   0.000   0.000
+    4ALA      N    5   0.600   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_awwa.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_awwa.gro
new file mode 100644
index 0000000000..7e7f6e5f92
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_awwa.gro
@@ -0,0 +1,10 @@
+Test gro file
+7
+    1ALA      O    1   0.000   0.000   0.000
+    2SOL     OW    2   0.300   0.000   0.000
+    2SOL    HW1    3   0.200   0.000   0.000
+    2SOL    HW2    4   0.400   0.000   0.000
+    3SOL     OW    5   0.600   0.000   0.000
+    3SOL    HW1    6   0.700   0.000   0.000
+    4ALA      O    7   0.900   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_awwwa.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_awwwa.gro
new file mode 100644
index 0000000000..78a065cf9d
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_awwwa.gro
@@ -0,0 +1,12 @@
+Test gro file
+9
+    1ALA      O    1   0.000   0.000   0.000
+    2SOL     OW    2   0.300   0.000   0.000
+    2SOL    HW1    3   0.200   0.000   0.000
+    2SOL    HW2    4   0.400   0.000   0.000
+    3SOL     OW    5   0.600   0.000   0.000
+    3SOL    HW1    6   0.700   0.000   0.000
+    4SOL     OW    7   0.900   0.000   0.000
+    4SOL    HW1    8   1.000   0.000   0.000
+    5ALA      O    9   1.200   0.000   0.000
+  10.0   10.0   10.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_awwwwa.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_awwwwa.gro
new file mode 100644
index 0000000000..22068b5f7b
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_awwwwa.gro
@@ -0,0 +1,14 @@
+Test gro file
+11
+    1ALA      O    1   0.000   0.000   0.000
+    2SOL     OW    2   0.300   0.000   0.000
+    2SOL    HW1    3   0.200   0.000   0.000
+    2SOL    HW2    4   0.400   0.000   0.000
+    3SOL     OW    5   0.600   0.000   0.000
+    3SOL    HW1    6   0.700   0.000   0.000
+    4SOL     OW    7   0.900   0.000   0.000
+    4SOL    HW1    8   1.000   0.000   0.000
+    5SOL     OW    9   1.200   0.000   0.000
+    5SOL    HW1   10   1.300   0.000   0.000
+    6ALA      O   11   1.400   0.000   0.000
+  10.0   10.0   10.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_branch.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_branch.gro
new file mode 100644
index 0000000000..a488330fb3
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_branch.gro
@@ -0,0 +1,12 @@
+Test gro file
+9
+    1ALA      O    1   0.000   0.000   0.000
+    2SOL     OW    2   0.300   0.000   0.000
+    2SOL    HW1    3   0.200   0.000   0.000
+    2SOL    HW2    4   0.400   0.000   0.000
+    3SOL     OW    5   0.600   0.000   0.000
+    3SOL    HW1    6   0.700   0.000   0.000
+    3SOL    HW2    7   0.600   0.100   0.000
+    4ALA      O    8   0.900   0.000   0.000
+    5ALA      O    9   0.600   0.300   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_da.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_da.gro
new file mode 100644
index 0000000000..4ba4307698
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_da.gro
@@ -0,0 +1,6 @@
+Test gro file
+3
+    1ALA      N    1   0.000   0.000   0.000
+    1ALA      H    2   0.100   0.000   0.000
+    4ALA      O    3   0.300   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_da_pbc.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_da_pbc.gro
new file mode 100644
index 0000000000..f231b93abb
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_da_pbc.gro
@@ -0,0 +1,6 @@
+Test gro file
+3
+    1ALA      N    1   0.800   0.000   0.000
+    1ALA      H    2   0.900   0.000   0.000
+    4ALA      O    3   0.100   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_duplicate_water.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_duplicate_water.gro
new file mode 100644
index 0000000000..daf237d7d1
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_duplicate_water.gro
@@ -0,0 +1,10 @@
+Test gro file
+    7
+    1LEU      O    1   1.876   0.810   1.354
+  117SOL    HW1    2   1.853   0.831   1.162
+  117SOL     OW    3   1.877   0.890   1.081
+  117SOL    HW2    4   1.908   0.828   1.007
+  135SOL     OW    5   1.924   0.713   0.845
+    1LEU      H    6   1.997   0.991   1.194
+    1LEU      N    7   2.041   1.030   1.274
+   2.22092   2.22092   2.22092
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_dwa.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_dwa.gro
new file mode 100644
index 0000000000..76849bce17
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_dwa.gro
@@ -0,0 +1,8 @@
+Test gro file
+5
+    1ALA      N    1   0.000   0.000   0.000
+    1ALA      H    2   0.100   0.000   0.000
+    2SOL     OW    3   0.300   0.000   0.000
+    2SOL    HW2    4   0.400   0.000   0.000
+    4ALA      O    5   0.600   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_dwd.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_dwd.gro
new file mode 100644
index 0000000000..b0f864e488
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_dwd.gro
@@ -0,0 +1,8 @@
+Test gro file
+5
+    1ALA      N    1   0.000   0.000   0.000
+    1ALA      H    2   0.100   0.000   0.000
+    2SOL     OW    3   0.300   0.000   0.000
+    4ALA      H    4   0.500   0.000   0.000
+    4ALA      N    5   0.600   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_empty.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_empty.gro
new file mode 100644
index 0000000000..d920f81fc2
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_empty.gro
@@ -0,0 +1,8 @@
+Test gro file
+5
+    1ALA      N    1   0.000   0.000   0.000
+    1ALA      H    2   0.100   0.000   0.000
+    2SOL     OW    3   3.000   0.000   0.000
+    4ALA      H    4   0.500   0.000   0.000
+    4ALA      N    5   0.600   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_loop.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_loop.gro
new file mode 100644
index 0000000000..2cbb57eb74
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_loop.gro
@@ -0,0 +1,8 @@
+Test gro file
+5
+    1ALA      O    1   0.000   0.001   0.000
+    2SOL     OW    2   0.300   0.001   0.000
+    2SOL    HW1    3   0.200   0.002   0.000
+    2SOL    HW2    4   0.200   0.000   0.000
+    4ALA      O    5   0.600   0.000   0.000
+  1.0   1.0   1.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_multiframe.dcd b/testsuite/MDAnalysisTests/data/waterbridge/wb_multiframe.dcd
new file mode 100644
index 0000000000..88c6cf130a
Binary files /dev/null and b/testsuite/MDAnalysisTests/data/waterbridge/wb_multiframe.dcd differ
diff --git a/testsuite/MDAnalysisTests/data/waterbridge/wb_multiframe.gro b/testsuite/MDAnalysisTests/data/waterbridge/wb_multiframe.gro
new file mode 100644
index 0000000000..e69670fd60
--- /dev/null
+++ b/testsuite/MDAnalysisTests/data/waterbridge/wb_multiframe.gro
@@ -0,0 +1,16 @@
+Test gro file
+ 13
+    1ALA      O    1   0.000   0.000   0.000
+    1ALA      H    2   0.000   0.000   0.000
+    2SOL     OW    3   0.300   0.000   0.000
+    2SOL    HW1    4   0.200   0.000   0.000
+    2SOL    HW2    5   0.400   0.000   0.000
+    3SOL     OW    6   0.600   0.000   0.000
+    3SOL    HW1    7   0.700   0.000   0.000
+    4SOL     OW    8   0.900   0.000   0.000
+    4SOL    HW1    9   1.000   0.000   0.000
+    5SOL     OW   10   1.200   0.000   0.000
+    5SOL    HW1   11   1.300   0.000   0.000
+    6ALA      H   12   1.400   0.000   0.000
+    6ALA      O   13   1.400   0.000   0.000
+   10.0   10.0   10.0
\ No newline at end of file
diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py
index 396bdbd86a..0c2f085d7b 100644
--- a/testsuite/MDAnalysisTests/datafiles.py
+++ b/testsuite/MDAnalysisTests/datafiles.py
@@ -396,8 +396,24 @@
     "SURFACE_PDB",  # 111 FCC lattice topology for NSGrid bug #2345
     "SURFACE_TRR",  # full precision coordinates for NSGrid bug #2345
     "DSSP",  # DSSP test suite
+    "WB_AD",  # GRO file with one hydrogen bond donor bonding to a hydrogen bond acceptor
+    "WB_AWA",  # GRO file with two hydrogen bond acceptors that are joined by a water
+    "WB_AWA_AWWA",  # GRO file with one hydrogen bond acceptors are bonded through one or two water
+    "WB_AWD",  # GRO file with one hydrogen bond acceptor bonding to a hydrogen bond donor through a water
+    "WB_AWWA",  # GRO file with one hydrogen bond acceptor bonding to a hydrogen bond acceptor through two waters
+    "WB_AWWWA",  # GRO file with one hydrogen bond acceptor bonding to a hydrogen bond acceptor through three waters
+    "WB_AWWWWA",  # GRO file with one hydrogen bond acceptor bonding to a hydrogen bond acceptor through four waters
+    "WB_BRANCH",  # GRO file with one hydrogen bond acceptor bonding to two hydrogen bond acceptor in selection 2
+    "WB_DA",  # GRO file with one hydrogen bond acceptor bonding to a hydrogen bond donor
+    "WB_DA_PBC",  # GRO file with one hydrogen bond acceptor bonding to a hydrogen bond donor but in a PBC condition
+    "WB_DWA",  # GRO file with one hydrogen bond donor bonding to a hydrogen bond acceptor through a water
+    "WB_DWD",  # GRO file with one hydrogen bond donor bonding to a hydrogen bond donor through a water
+    "WB_EMPTY",  # GRO file with no hydrogen bonds
+    "WB_LOOP",  # GRO with hydrogen bond acceptor bonding to water which bonds back to the first hydrogen bond acceptor to form a loop
+    "WB_DUPLICATE_WATER",  # GRO file with duplicate waters, reference to case #3119
+    "WB_MULTIFRAME_GRO",  # topology file of containing water bridges for multiple frame tests
+    "WB_MULTIFRAME_DCD",  # trajectory file of containing water bridges for multiple frame tests
 ]
-
 from importlib import resources
 import MDAnalysisTests.data
 
@@ -928,5 +944,27 @@
 # DSSP testing: from https://github.com/ShintaroMinami/PyDSSP
 DSSP = (_data_ref / "dssp").as_posix()
 
+
+# Water Bridge files. Currently used due to the
+# incompatibility of StringIO with parallelization backends. Issue #5221
+WB_AD = (_data_ref / "waterbridge/wb_ad.gro").as_posix()
+WB_AWA = (_data_ref / "waterbridge/wb_awa.gro").as_posix()
+WB_AWA_AWWA = (_data_ref / "waterbridge/wb_awa_awwa.gro").as_posix()
+WB_AWD = (_data_ref / "waterbridge/wb_awd.gro").as_posix()
+WB_AWWA = (_data_ref / "waterbridge/wb_awwa.gro").as_posix()
+WB_AWWWA = (_data_ref / "waterbridge/wb_awwwa.gro").as_posix()
+WB_AWWWWA = (_data_ref / "waterbridge/wb_awwwwa.gro").as_posix()
+WB_BRANCH = (_data_ref / "waterbridge/wb_branch.gro").as_posix()
+WB_DA = (_data_ref / "waterbridge/wb_da.gro").as_posix()
+WB_DA_PBC = (_data_ref / "waterbridge/wb_da_pbc.gro").as_posix()
+WB_DWA = (_data_ref / "waterbridge/wb_dwa.gro").as_posix()
+WB_DWD = (_data_ref / "waterbridge/wb_dwd.gro").as_posix()
+WB_EMPTY = (_data_ref / "waterbridge/wb_empty.gro").as_posix()
+WB_LOOP = (_data_ref / "waterbridge/wb_loop.gro").as_posix()
+WB_DUPLICATE_WATER = (
+    _data_ref / "waterbridge/wb_duplicate_water.gro"
+).as_posix()
+WB_MULTIFRAME_GRO = (_data_ref / "waterbridge/wb_multiframe.gro").as_posix()
+WB_MULTIFRAME_DCD = (_data_ref / "waterbridge/wb_multiframe.dcd").as_posix()
 # This should be the last line: clean up namespace
 del resources