diff --git a/.github/workflows/full_test.yml b/.github/workflows/full_test.yml
index 26a261d3a..256d06d7e 100644
--- a/.github/workflows/full_test.yml
+++ b/.github/workflows/full_test.yml
@@ -17,7 +17,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, windows-latest]
-        python-version: [3.8, "3.12"]
+        python-version: [3.9, "3.12"]
 
     steps:
     - uses: actions/checkout@v3
diff --git a/.github/workflows/legacy_test.yml b/.github/workflows/legacy_test.yml
index 63f746dce..77ac50f87 100644
--- a/.github/workflows/legacy_test.yml
+++ b/.github/workflows/legacy_test.yml
@@ -17,7 +17,7 @@ jobs:
         - name: Set up Python
           uses: actions/setup-python@v4
           with:
-            python-version: 3.8
+            python-version: 3.9
 
         - name: Install older dependencies
           run: |
diff --git a/mikeio/dataset/_dataarray.py b/mikeio/dataset/_dataarray.py
index 556f1552e..f32425b93 100644
--- a/mikeio/dataset/_dataarray.py
+++ b/mikeio/dataset/_dataarray.py
@@ -1873,6 +1873,54 @@ def to_xarray(self) -> "xarray.DataArray":
         )
         return xr_da
 
+    def to_xugrid(self):
+        import xugrid as xu
+        import xarray as xr
+
+        assert isinstance(
+            self.geometry, GeometryFM2D
+        ), "Only flexible mesh 2D data is supported"
+        # face_code_connectivity is slightly different than the element_table
+
+        g = self.geometry
+
+        if g.is_tri_only:
+            conn = np.zeros(shape=(self.geometry.n_elements, 3), dtype=np.int32)
+            for i in range(g.n_elements):
+                conn[i, 0:3] = g.element_table[i]
+        else:
+            conn = np.zeros(shape=(self.geometry.n_elements, 4), dtype=np.int32)
+            conn[:] = -1
+
+            for i in range(g.n_elements):
+                if len(g.element_table[i]) == 3:
+                    conn[i, 0:3] = g.element_table[i]
+                else:
+                    conn[i] = g.element_table[i]
+
+        xn = g.node_coordinates[:, 0]
+        yn = g.node_coordinates[:, 1]
+
+        grid = xu.Ugrid2d(xn, yn, fill_value=-1, face_node_connectivity=conn)
+
+        if "time" not in self.dims:
+            da = xr.DataArray(
+                name=self.name,
+                data=self.to_numpy(),
+                coords={},
+                dims=[grid.face_dimension],
+            )
+        else:
+            da = xr.DataArray(
+                name=self.name,
+                data=self.to_numpy(),
+                coords={"time": self.time},
+                dims=["time", grid.face_dimension],
+            )
+
+        uda = xu.UgridDataArray(da, grid)
+        return uda
+
     # ===============================================
 
     def __repr__(self) -> str:
diff --git a/mikeio/dataset/_dataset.py b/mikeio/dataset/_dataset.py
index 00e1b1974..b823fb998 100644
--- a/mikeio/dataset/_dataset.py
+++ b/mikeio/dataset/_dataset.py
@@ -1898,6 +1898,15 @@ def to_xarray(self) -> "xarray.Dataset":
         data = {da.name: da.to_xarray() for da in self}
         return xarray.Dataset(data)
 
+    def to_xugrid(self):
+        import xugrid as xu
+
+        data = [da.to_xugrid() for da in self]
+        ds = xu.UgridDataset(grids=data[0].grid)
+        for da in data:
+            ds[da.name] = da
+        return ds
+
     # ===============================================
 
     def __repr__(self) -> str:
diff --git a/mypy.ini b/mypy.ini
index 385b74d70..74d5feceb 100644
--- a/mypy.ini
+++ b/mypy.ini
@@ -33,4 +33,7 @@ ignore_missing_imports = True
 [mypy-yaml.*]
 ignore_missing_imports = True
 
+[mypy-xugrid.*]
+ignore_missing_imports = True
+
 
diff --git a/pyproject.toml b/pyproject.toml
index 795bf08c7..86886c0bb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -58,7 +58,7 @@ dev = ["pytest",
        "mypy==1.6.1",
        ]
 
-test = ["pytest", "pytest-cov", "xarray","mypy==1.6.1","shapely","pyproj"]
+test = ["pytest", "pytest-cov", "xarray","mypy==1.6.1","shapely","pyproj","xugrid"]
 
 notebooks= [
             "nbformat",
diff --git a/tests/test_dfsu.py b/tests/test_dfsu.py
index f3092f907..154ee153c 100644
--- a/tests/test_dfsu.py
+++ b/tests/test_dfsu.py
@@ -970,3 +970,29 @@ def test_writing_non_equdistant_dfsu_is_not_possible(tmp_path):
     with pytest.raises(ValueError, match="equidistant"):
         fp = tmp_path / "not_gonna_work.dfsu"
         dss.to_dfs(fp)
+
+
+def test_convert_dfsu2d_to_xugrid(tmp_path):
+
+    import xugrid
+
+    ds = mikeio.read("tests/testdata/NorthSea_HD_and_windspeed.dfsu")
+    ds.sel(y=55, x=0).isel(time=-1).to_numpy() == pytest.approx(0.25443718)
+    xu_ds = ds.to_xugrid()
+
+    xu_ds.ugrid.sel(y=55, x=0)["Surface elevation"].isel(
+        time=-1
+    ).to_numpy() == pytest.approx(0.25443718)
+    assert xu_ds.time[-1] == ds.time[-1]
+
+    fp = tmp_path / "north_sea.nc"
+    xu_ds.ugrid.to_netcdf(fp)
+
+    xu_ds2 = xugrid.open_dataset(fp)
+    assert xu_ds2.time[-1] == ds.time[-1]
+
+
+def test_mixed_mesh_time_invariant_to_xugrid():
+    ds = mikeio.read("tests/testdata/FakeLake.dfsu", time=0)
+    xu_ds = ds.to_xugrid()
+    assert len(xu_ds.mesh2d_nFaces) == 1011