diff --git a/README.md b/README.md
index d3b8fba3..02ce6353 100644
--- a/README.md
+++ b/README.md
@@ -125,7 +125,9 @@ QGIS contains internal copy of MDAL library in following versions:
 | 3.30.0  | 1.0.2   | |
 | 3.36.0  | 1.1.0   | Mike21 format support read/write |
 | 3.38.0  | 1.2.0   | Groundwater / surface water meshes for 3Di format |
-| 3.42.0  | 1.3.0   | Fix 2dm format coordinates saving, Support Dataset Group Removal From Mesh |
+| 3.42.0  | 1.3.0   | Fix 2dm format coordinates saving |
+|         |         | XMDF loading Dataset Group fix, support for Mesh in XMDF files (as 2DMeshModule) |
+|         |         | Support Dataset Group Removal From Mesh |
 
 versions `X.Y.9Z` are development versions or alpha/beta releases (e.g. `0.4.90`, `0.4.91`, ...)
 
diff --git a/docs/source/drivers/xmdf.rst b/docs/source/drivers/xmdf.rst
index 851bc198..da964089 100644
--- a/docs/source/drivers/xmdf.rst
+++ b/docs/source/drivers/xmdf.rst
@@ -10,5 +10,8 @@ XMDF -- eXtensible Model Data Format
 
 MDAL supports reading of the XMDF format generated by TUFLOW_, HYDRO_AS-2D_ and other hydraulic modelling software applications.
 
+Since 1.3 MDAL supports reading mesh stored in XMDF format according to specification - http://xmdf.aquaveo.com/doc2.2/html/modules.html. 
+
 .. _TUFLOW: https://www.tuflow.com/
 .. _HYDRO_AS-2D : https://www.hydroas-2d.com/
+.. _Mesh Format: http://xmdf.aquaveo.com/doc2.2/html/group__d2d4d1d.html
\ No newline at end of file
diff --git a/mdal/frmts/mdal_xmdf.cpp b/mdal/frmts/mdal_xmdf.cpp
index c09d5fa0..340722b9 100644
--- a/mdal/frmts/mdal_xmdf.cpp
+++ b/mdal/frmts/mdal_xmdf.cpp
@@ -92,8 +92,8 @@ size_t MDAL::XmdfDataset::activeData( size_t indexStart, size_t count, int *buff
 MDAL::DriverXmdf::DriverXmdf()
   : Driver( "XMDF",
             "TUFLOW XMDF",
-            "*.xmdf",
-            Capability::ReadDatasets )
+            "*.xmdf;;*.h5",
+            Capability::ReadDatasets | Capability::ReadMesh )
 {
 }
 
@@ -206,6 +206,17 @@ void MDAL::DriverXmdf::addDatasetGroupsFromXmdfGroup( DatasetGroups &groups,
     size_t vertexCount,
     size_t faceCount ) const
 {
+  // check if this root group can be loaded as a dataset group and if so, then load it
+  std::vector<std::string> gDataNames = rootGroup.datasets();
+  if ( MDAL::contains( gDataNames, "Times" ) &&
+       MDAL::contains( gDataNames, "Values" ) &&
+       MDAL::contains( gDataNames, "Mins" ) &&
+       MDAL::contains( gDataNames, "Maxs" ) )
+  {
+    std::shared_ptr<DatasetGroup> ds = readXmdfGroupAsDatasetGroup( rootGroup, rootGroup.name() + nameSuffix, vertexCount, faceCount );
+    groups.push_back( ds );
+  }
+
   for ( const std::string &groupName : rootGroup.groups() )
   {
     HdfGroup g = rootGroup.group( groupName );
@@ -331,3 +342,241 @@ std::shared_ptr<MDAL::DatasetGroup> MDAL::DriverXmdf::readXmdfGroupAsDatasetGrou
 
   return group;
 }
+
+bool MDAL::DriverXmdf::canReadMesh( const std::string &uri )
+{
+  HdfFile file( uri, HdfFile::ReadOnly );
+  if ( !file.isValid() )
+  {
+    return false;
+  }
+
+  HdfDataset dsFileType = file.dataset( "/File Type" );
+  if ( dsFileType.readString() != "Xmdf" )
+  {
+    return false;
+  }
+
+  std::vector<std::string> meshPaths = meshGroupPaths( file );
+
+  return !meshPaths.empty();
+}
+
+std::string MDAL::DriverXmdf::buildUri( const std::string &meshFile )
+{
+  mDatFile = meshFile;
+
+  std::vector<std::string> meshNames = findMeshesNames();
+
+  return MDAL::buildAndMergeMeshUris( meshFile, meshNames, name() );
+}
+
+std::vector<std::string> MDAL::DriverXmdf::findMeshesNames() const
+{
+  std::vector<std::string> meshesInFile;
+
+  HdfFile file( mDatFile, HdfFile::ReadOnly );
+  if ( !file.isValid() )
+  {
+    return meshesInFile;
+  }
+
+  meshesInFile = meshGroupPaths( file );
+
+  return meshesInFile;
+}
+
+std::vector<std::string> MDAL::DriverXmdf::meshGroupPaths( const HdfFile &file ) const
+{
+  std::vector<std::string> meshPaths;
+
+  std::vector<std::string> rootGroups = file.groups();
+
+  for ( const std::string &groupName : rootGroups )
+  {
+    HdfGroup g = file.group( groupName );
+    std::vector<std::string> paths = meshGroupPaths( g );
+    meshPaths.insert( meshPaths.end(), paths.begin(), paths.end() );
+  }
+
+  return meshPaths;
+}
+
+std::vector<std::string> MDAL::DriverXmdf::meshGroupPaths( const HdfGroup &group ) const
+{
+  std::vector<std::string> meshPaths;
+
+  std::vector<std::string> gDataNames = group.groups();
+
+  if ( MDAL::contains( gDataNames, "Nodes" ) ||
+       MDAL::contains( gDataNames, "Elements" ) )
+  {
+    meshPaths.push_back( group.name() );
+  }
+
+  for ( const std::string &groupName : gDataNames )
+  {
+    HdfGroup g = group.group( groupName );
+    std::vector<std::string> paths = meshGroupPaths( g );
+    meshPaths.insert( meshPaths.end(), paths.begin(), paths.end() );
+  }
+
+  return meshPaths;
+}
+
+std::unique_ptr< MDAL::Mesh > MDAL::DriverXmdf::load( const std::string &meshFile, const std::string &meshName )
+{
+  mDatFile = meshFile;
+
+  MDAL::Log::resetLastStatus();
+
+  HdfFile file( mDatFile, HdfFile::ReadOnly );
+  if ( !file.isValid() )
+  {
+    MDAL::Log::error( MDAL_Status::Err_UnknownFormat, name(), "File " + mDatFile + " is not valid" );
+    return nullptr;
+  }
+
+  HdfDataset dsFileType = file.dataset( "/File Type" );
+  if ( dsFileType.readString() != "Xmdf" )
+  {
+    MDAL::Log::error( MDAL_Status::Err_UnknownFormat, name(), "Unknown dataset file type" );
+    return nullptr;
+  }
+
+  std::vector<std::string> meshNames = findMeshesNames();
+
+  if ( meshNames.empty() )
+  {
+    MDAL::Log::error( MDAL_Status::Err_IncompatibleMesh, name(), "No meshes found in file " + mDatFile );
+    return nullptr;
+  }
+
+  std::string meshNameToLoad = meshName;
+
+  if ( meshNameToLoad.empty() )
+  {
+    meshNameToLoad = meshNames[0];
+  }
+
+  if ( !MDAL::contains( meshNames, meshNameToLoad ) )
+  {
+    MDAL::Log::error( MDAL_Status::Err_IncompatibleMesh, name(), "No meshes with name " + meshNameToLoad + " found in file " + mDatFile );
+    return nullptr;
+  }
+
+  HdfGroup groupMeshModule = file.group( meshNameToLoad );
+
+  std::vector<std::string> gDataNames = groupMeshModule.groups();
+
+  HdfGroup gNodes = groupMeshModule.group( "Nodes" );
+
+  std::vector<std::string> namesNodes = gNodes.datasets();
+  HdfDataset nodes = gNodes.dataset( namesNodes[0] );
+
+  std::vector<hsize_t> nodesDims = nodes.dims();
+  hsize_t nodesRows = nodesDims[0];
+  size_t vertexDims = nodesDims[1];
+
+  if ( vertexDims < 2 || vertexDims > 3 )
+  {
+    MDAL::Log::error( MDAL_Status::Err_IncompatibleMesh, name(), "Vertices have unsupported number of dimensions " + std::to_string( vertexDims ) + " only 2 (X,Y) or 3 (X, Y, Z) dimensions are allowed." );
+    return nullptr;
+  }
+
+  std::vector<double> nodesData = nodes.readArrayDouble();
+
+  Vertices vertices( nodesRows );
+
+  size_t currentVertexIndex = 0;
+  size_t i = 0;
+  while ( i < nodesData.size() )
+  {
+    Vertex &vertex = vertices[currentVertexIndex];
+
+    vertex.x = nodesData[i];
+    i++;
+    vertex.y = nodesData[i];
+    i++;
+    if ( vertexDims == 3 )
+    {
+      vertex.z = nodesData[i];
+      i++;
+    }
+    currentVertexIndex++;
+  }
+
+  nodesData.clear();
+
+  HdfGroup gElements = groupMeshModule.group( "Elements" );
+
+  std::vector<std::string> namesElements = gElements.datasets();
+  HdfDataset elements = gElements.dataset( namesElements[0] );
+
+  std::vector<hsize_t> elementsDims = elements.dims();
+  hsize_t elementsRows = elementsDims[0];
+  int elementsRowsDims = elementsDims[1];
+
+  std::vector<int> facesData = elements.readArrayInt();
+
+  Faces faces( elementsRows );
+  int maxVerticesPerFace = 0;
+
+  size_t currentFaceIndex = 0;
+  i = 0;
+  while ( i < facesData.size() )
+  {
+    std::vector<size_t> tempFace;
+    for ( int j = 0; j < elementsRowsDims; j++ )
+    {
+      int vertexIndex = facesData[i];
+      if ( vertexIndex > 0 )
+      {
+        // XMDF is 1-based, MDAL is 0-based
+        tempFace.push_back( facesData[i] - 1 );
+      }
+      i++;
+    }
+
+    // only store faces with more than 2 vertices
+    if ( tempFace.size() > 2 )
+    {
+      Face &face = faces[currentFaceIndex];
+      std::copy( tempFace.begin(), tempFace.end(), std::back_inserter( face ) );
+
+      if ( tempFace.size() > maxVerticesPerFace )
+      {
+        maxVerticesPerFace = tempFace.size();
+      }
+
+      currentFaceIndex++;
+    }
+  }
+
+  facesData.clear();
+
+  // copy only the faces that have been properly filled
+  faces = Faces( faces.begin(), faces.begin() + currentFaceIndex );
+
+  // create the mesh and set the required data
+  std::unique_ptr< MemoryMesh > mesh(
+    new MemoryMesh(
+      name(),
+      maxVerticesPerFace,
+      mDatFile
+    )
+  );
+
+  std::vector<double> values( vertices.size() );
+  for ( size_t i = 0; i < vertices.size(); ++i )
+  {
+    values[i] = vertices[i].z;
+  }
+
+  mesh->setFaces( std::move( faces ) );
+  mesh->setVertices( std::move( vertices ) );
+
+  addVertexScalarDatasetGroup( mesh.get(), values, "Z-Values" );
+
+  return mesh;
+}
diff --git a/mdal/frmts/mdal_xmdf.hpp b/mdal/frmts/mdal_xmdf.hpp
index cd565034..3163031b 100644
--- a/mdal/frmts/mdal_xmdf.hpp
+++ b/mdal/frmts/mdal_xmdf.hpp
@@ -89,6 +89,9 @@ namespace MDAL
       bool canReadDatasets( const std::string &uri ) override;
       void load( const std::string &datFile, Mesh *mesh ) override;
 
+      bool canReadMesh( const std::string &uri ) override;
+      std::unique_ptr< Mesh > load( const std::string &meshFile, const std::string &meshName = "" ) override;
+
     private:
       MDAL::Mesh *mMesh = nullptr;
       std::string mDatFile;
@@ -111,6 +114,11 @@ namespace MDAL
                            size_t vertexCount,
                            size_t faceCount ) const;
 
+      std::string buildUri( const std::string &meshFile ) override;
+      std::vector<std::string> findMeshesNames() const;
+
+      std::vector<std::string> meshGroupPaths( const HdfGroup &group ) const;
+      std::vector<std::string> meshGroupPaths( const HdfFile &file ) const;
   };
 
 } // namespace MDAL
diff --git a/tests/data/xmdf/withMesh/data.h5 b/tests/data/xmdf/withMesh/data.h5
new file mode 100644
index 00000000..51ebdc14
Binary files /dev/null and b/tests/data/xmdf/withMesh/data.h5 differ
diff --git a/tests/data/xmdf/withMesh/mesh.h5 b/tests/data/xmdf/withMesh/mesh.h5
new file mode 100644
index 00000000..48a7ad5c
Binary files /dev/null and b/tests/data/xmdf/withMesh/mesh.h5 differ
diff --git a/tests/data/xmdf/withMesh/multiple_meshes.h5 b/tests/data/xmdf/withMesh/multiple_meshes.h5
new file mode 100644
index 00000000..757669ab
Binary files /dev/null and b/tests/data/xmdf/withMesh/multiple_meshes.h5 differ
diff --git a/tests/test_xmdf.cpp b/tests/test_xmdf.cpp
index c275a259..18de75ec 100644
--- a/tests/test_xmdf.cpp
+++ b/tests/test_xmdf.cpp
@@ -501,6 +501,150 @@ TEST( MeshXmdfTest, withFinalgroup )
   MDAL_CloseMesh( m );
 }
 
+TEST( MeshXmdfTest, MeshLoading )
+{
+  std::string pathMesh = test_file( "/xmdf/withMesh/mesh.h5" );
+  MDAL_MeshH m = MDAL_LoadMesh( pathMesh.c_str() );
+  ASSERT_NE( m, nullptr );
+
+  ASSERT_EQ( 5, MDAL_M_vertexCount( m ) );
+  ASSERT_EQ( 2, MDAL_M_faceCount( m ) );
+
+  double minX, maxX, minY, maxY;
+  MDAL_M_extent( m, &minX, &maxX, &minY, &maxY );
+
+  EXPECT_DOUBLE_EQ( 1000, minX );
+  EXPECT_DOUBLE_EQ( 3000, maxX );
+  EXPECT_DOUBLE_EQ( 2000, minY );
+  EXPECT_DOUBLE_EQ( 3000, maxY );
+
+  EXPECT_EQ( 1, MDAL_M_datasetGroupCount( m ) );
+
+  MDAL_DatasetGroupH g = MDAL_M_datasetGroup( m, 0 );
+  ASSERT_NE( g, nullptr );
+
+  EXPECT_TRUE( MDAL_G_hasScalarData( g ) );
+
+  int meta_count = MDAL_G_metadataCount( g );
+  ASSERT_EQ( 1, meta_count );
+
+  const char *name = MDAL_G_name( g );
+  EXPECT_EQ( std::string( "Z-Values" ), std::string( name ) );
+
+  EXPECT_EQ( MDAL_M_datasetGroupCount( m ), 1 );
+
+  MDAL_DatasetH ds = MDAL_G_dataset( g, 0 );
+
+  EXPECT_TRUE( MDAL_D_isValid( ds ) );
+
+  double min, max;
+  MDAL_D_minimumMaximum( ds, &min, &max );
+
+  EXPECT_DOUBLE_EQ( 10, min );
+  EXPECT_DOUBLE_EQ( 50, max );
+}
+
+TEST( MeshXmdfTest, DataInSubfolder )
+{
+  std::string pathMesh = test_file( "/xmdf/withMesh/mesh.h5" );
+  MDAL_MeshH m = MDAL_LoadMesh( pathMesh.c_str() );
+  ASSERT_NE( m, nullptr );
+  ASSERT_EQ( 1, MDAL_M_datasetGroupCount( m ) );
+
+  std::string path = test_file( "/xmdf/withMesh/data.h5" );
+  EXPECT_TRUE( std::string( MDAL_MeshNames( path.c_str() ) ).empty() );
+  MDAL_M_LoadDatasets( m, path.c_str() );
+  MDAL_Status s = MDAL_LastStatus();
+  EXPECT_EQ( MDAL_Status::None, s );
+
+  ASSERT_EQ( 3, MDAL_M_datasetGroupCount( m ) );
+
+  MDAL_DatasetGroupH g = MDAL_M_datasetGroup( m, 1 );
+  ASSERT_NE( g, nullptr );
+  ASSERT_EQ( std::string( "/Datasets/data" ), std::string( MDAL_G_name( g ) ) );
+
+  MDAL_DatasetH ds = MDAL_G_dataset( g, 0 );
+  ASSERT_NE( ds, nullptr );
+  ASSERT_TRUE( MDAL_D_isValid( ds ) );
+  EXPECT_EQ( 5, MDAL_D_valueCount( ds ) );
+
+  double min, max;
+  MDAL_D_minimumMaximum( ds, &min, &max );
+  EXPECT_DOUBLE_EQ( 1, min );
+  EXPECT_DOUBLE_EQ( 3, max );
+
+  g = MDAL_M_datasetGroup( m, 2 );
+  ASSERT_NE( g, nullptr );
+  ASSERT_EQ( std::string( "subgroup" ), std::string( MDAL_G_name( g ) ) );
+
+  ds = MDAL_G_dataset( g, 0 );
+  ASSERT_NE( ds, nullptr );
+  ASSERT_TRUE( MDAL_D_isValid( ds ) );
+  EXPECT_EQ( 5, MDAL_D_valueCount( ds ) );
+
+  MDAL_D_minimumMaximum( ds, &min, &max );
+  EXPECT_DOUBLE_EQ( 100, min );
+  EXPECT_DOUBLE_EQ( 500, max );
+}
+
+TEST( MeshXmdfTest, MultipleMeshes )
+{
+  // dataset without mesh
+  std::string path = test_file( "/xmdf/withMesh/data.h5" );
+  EXPECT_EQ( std::string( MDAL_MeshNames( path.c_str() ) ), std::string( "" ) );
+  MDAL_MeshH m = MDAL_LoadMesh( path.c_str() );
+  ASSERT_EQ( m, nullptr );
+
+  // test that all 4 meshes are found
+  path = test_file( "/xmdf/withMesh/multiple_meshes.h5" );
+  EXPECT_EQ( MDAL_MeshNames( path.c_str() ),
+             "XMDF:\"" + path + "\":/2DMeshModule/triangle_and_quad;;" +
+             "XMDF:\"" + path + "\":/AnotherGroup/2DMeshModule/triangle_and_quad;;" +
+             "XMDF:\"" + path + "\":/OtherModuleWithMesh/t_q;;" +
+             "XMDF:\"" + path + "\":/YetAnotherGroup/ModuleWithMesh/t_q" );
+
+  // default mesh can be loaded
+  m = MDAL_LoadMesh( path.c_str() );
+  ASSERT_NE( m, nullptr );
+
+  // load specific mesh and test extent to validate correct mesh was loaded
+  std::string uri = "\"" + path + "\":/2DMeshModule/triangle_and_quad";
+  m = MDAL_LoadMesh( uri.c_str() );
+  ASSERT_NE( m, nullptr );
+
+  double minX, maxX, minY, maxY;
+  MDAL_M_extent( m, &minX, &maxX, &minY, &maxY );
+
+  EXPECT_DOUBLE_EQ( 1000, minX );
+  EXPECT_DOUBLE_EQ( 3000, maxX );
+  EXPECT_DOUBLE_EQ( 2000, minY );
+  EXPECT_DOUBLE_EQ( 3000, maxY );
+
+  // load another mesh
+  uri = "\"" + path + "\":/AnotherGroup/2DMeshModule/triangle_and_quad";
+  m = MDAL_LoadMesh( uri.c_str() );
+  ASSERT_NE( m, nullptr );
+
+  MDAL_M_extent( m, &minX, &maxX, &minY, &maxY );
+
+  EXPECT_DOUBLE_EQ( 10, minX );
+  EXPECT_DOUBLE_EQ( 30, maxX );
+  EXPECT_DOUBLE_EQ( 20, minY );
+  EXPECT_DOUBLE_EQ( 30, maxY );
+
+  // load another mesh
+  uri = "\"" + path + "\":/YetAnotherGroup/ModuleWithMesh/t_q";
+  m = MDAL_LoadMesh( uri.c_str() );
+  ASSERT_NE( m, nullptr );
+
+  MDAL_M_extent( m, &minX, &maxX, &minY, &maxY );
+
+  EXPECT_DOUBLE_EQ( 1, minX );
+  EXPECT_DOUBLE_EQ( 3, maxX );
+  EXPECT_DOUBLE_EQ( 2, minY );
+  EXPECT_DOUBLE_EQ( 3, maxY );
+}
+
 int main( int argc, char **argv )
 {
   testing::InitGoogleTest( &argc, argv );