I Charpentier
email: icharpentier@unistra.fr

XGeoTiles: Tiling watersheds with hexagons for serious games

Serious games on water governance at bring people on a particular question together such that they may overcome traditional positions and to trigger positive changes in both representation and behavior. To play with simplied but realistic watershed characteristics, we propose the XGeoTiles framework for the description of watershed elevation. It comes in the form of a 3D game board that is divided into hexagonal tiles (2D) and prisms (3D). Additional information such as land cover or stream network is also accounted.

Provided with a digital elevation model, land cover data and an outlet, the software returns a hexagonal model and a 3D OpenSCAD model that can be printed on appropriate devices. As an illustration, the hexagonal tiling is applied to a slightly undulating watershed. Results conrm a good representation of the hydrological and landscape properties when a relevant tile size is used.

Introduction

Serious games on water governance bring people on a particular question together such that they may overcome traditional positions and may trigger positive changes in both representation and behavior. Boards and role-playing, possibly enforced by computer simulations, have become key tools [START_REF] Medema | Understanding Game-based Approaches for Improving Sustainable Water Governance: The Potential of Serious Games to Solve Water Problems[END_REF] to study water management in the presence of stakeholders. The virtual river game [START_REF] Den Haan | The virtual river game: Gaming using models to collaboratively explore river management complexity[END_REF] is a good illustration of such a combination. Therein, a hexagonal tiling models rivers, channels and dikes. More generally, hexagonal boards, role play and simulation are also found in productivity and strategy game applications.

In the gaming framework, a hex tiling is a partition of a planar surface into hexagonal tiles. Present in nature-based solutions like honeycombs or crystals, hex tiles and hex prisms have a number of useful properties. Given a tile, the availability of 6 adjacent tiles and the equal distance between their centroids have retained the attention in a large number of strategy games. These tessellations oer more realistic movements on the board and a nicer representation of the landscape. Studied since a long time, the hexagonal tiling is receiving a renewed attention in the GIS modeling community. For instance, Wang et al. [START_REF] Wang | The isotropic organization of dem structure and extraction of valley lines using hexagonal grid[END_REF] model terrains from lidar data using a hexagonal lattice for the elevations to approximate the stream network along the edges.

Several methods are available to set up a 2D hex board. Hexagonal grid gaming papers sold for games orcid(s): 0000-0001-7511-2910 (I. Charpentier) such as Dungeons and Dragons may be used in a basic manner. More generally, software and apps (e.g. HeXMapMaker or Hextml) exist to create boards with landscape hexagonal tiles that describe cities, agriculture, forests, and water bodies. Some user-dened stream network may be added. To date, such boards remain user-dened 2D landscapes, with a reduced physical meaning. However, unsurprisingly, these key features are also present in the so-called Corine Land Cover (CLC) data sets [4].

The XGeoTiles framework is designed go one step beyond. Its goals are to account for watershed elevation and land cover in a simplied but realistic manner. A hexagonal tiling is implemented simultaneously for the digital elevation model (DEM) and the land cover data. Each of the hexagons is assimilated as a at horizontal homogeneous land parcel, the height and cover of which are the minimum height and the dominant cover over the parcel, respectively. This implementation yields a hex digital elevation model (XDEM) of the watershed made of hexagonal prisms. Classical methods for watershed delineation, ow direction, ow accumulation, drainage basin and stream computation are to be adapted to restore the hydrological functioning of the watershed. Furthermore, the tiling of land cover into hexagons (XLC) may provide information on some biophysical processes, fauna and ora habitats, and practices such as agriculture, urbanization and land management.

The simplied bio-geophysical information stored in both XDEM and XLC may be, rstly, processed to print 2D or 3D boards and, secondly, accounted into rainfall-runo simulations to play serious games with a realistic watershed board.

The layout of the paper is as follows. Some generalities about hexagon and square grids are recalled in section 2. The case study is introduced early in the paper for a better understanding. The workow and implementation details, including tile descriptors and data storage, are provided in section 3. Numerical results from dierent points of view are presented and discussed in Section 4. Finally, the main ndings of this study are summarized in the conclusion, section 5.

From terrain to landscape descriptors

A number of Digital Elevation Model (DEM) and land cover data (CLC, for instance) are nowadays freely available as rasters (tif images or ascii les) or shapeles, respectively. These may come with a dierent coordinate reference system (CRS) and/or a dierent resolution. A particular attention should be thus brought to their alignment, as discussed in subsection 2.1.

Some generalities about regular tilings are recalled in subsection 2.2 and some essential biogeophysical descriptors [START_REF] Grigal | Site and landscape characterization for ecological studies[END_REF] are introduced in subsection 2.3 prior to the design of boards from actual landscape datasets.

Case study and data preprocessing

The slightly undulating watershed of the Moderbach stream (Grand Est region, France) is chosen as an illustrative example [START_REF] Hariri | Pairing GIS and distributed hydrological models using Matlab 2[END_REF][START_REF] Viviani | Les étangs-réservoirs de la Ligne Maginot Aquatique, un socioécosystème durable ? 5ème colloque des Zones Ateliers-CNRS -2020[END_REF]. This stream has a length of 15 km and drains a rural catchment of 89 km 2 to the Albe stream, a left tributary of the Saar river. Elevations range from 210 to 305 m.a.s.l., see Fig. 1. The digital elevation model (DEM) raster is from BD Alti ®75m -V2. By denition, such a raster data model makes use of a uniform regular grid of square pixels. Raster are usually provided as tif images or ascii les. Land cover data provide tessellations/partitions into polygons that are either stored as shapeles or rasters. The land cover shapele used for Fig. 2 is from Corine Land Cover 12 with minimum mapping units of 0.25 km 2 . This rural landscape comprises a number of small villages (red shades, 10%), agricultural surfaces (yellow shades, 66%), forests (green shades, 22%), wetlands and 6 ponds (blue shades, 2 %) represented using polygons. Note that the rasterization of the polygons proposed in QGIS allows to create land cover grids. In the 1930s, the watershed was converted into a defense waterline based on 11 dams delineating 6 reservoirs (blue polygons of Fig. 2) on the Moderbach's tributaries and 5 forebays to be ooded along the main stream.

For consistency, DEM and land cover data has to be preprocessed to ensure the alignment of their CRS. Some reprojections and/or interpolations can be carried out on-the-y using tools such as ArcGIS or QGIS. Prospecting further watershed simulations, actual or spurious low areas of the DEM, the so-called sinks, that disturb the surface water ow should be preferably lled to ensure of a decrease of the heights from the most upstream cells to the outlet. Appendix A describes the data collection and the preprocessing steps for this case study. This may serve as a basis for other data sets and watersheds.

Regular tiling patterns

In the two-dimensional Euclidean plane, a tiling or tessellation of a surface consists of a partition into geometrical shapes, called tiles. Landscape datasets are good examples. On the one hand, see Fig. 1, the pixels of a raster image may be viewed as a regular tiling of the plane into small identical squares. On the other hand, the CLC polygons form an irregular tiling as in Fig. 2.

A regular tiling may be built by translating a reference square or a reference hexagon. Figure 3 presents square and hexagonal neighboring patterns, the tiles of which are of the same area. Both generation processes may be driven by the edge length and can be set up easily. The paper focuses on hexagonal tilings. The central reference hexagon, number 0, is dened by its centroid (x 0 , y 0 ) and the common apothem or inradius is hereafter denoted by a. The hexagon area is equal to

A = 2 √ 3a 2 . Its 6 vertices (x i 0 , y i 0 ) i=1,..,6
are equally distributed at a distance of e = 2

√ 3a/3
corresponding to the edge length, and an angle of 60°. Their coordinates satisfy:

x i = x 0 + a cos((2i -1)π/6), for i = 1, .., 6,
y i = y 0 + a sin((2i -1)π/6). (1)
The six neighboring tiles of the tile 0 are equally distributed using and α = 60°and dene a neighborhood with an area of 7A. In the local coordinate system centered on (x 0 , y 0 ), neighbors centroids are located at points:

x j = x 0 + 2a cos(jπ/3), for j = 1, .., 6, y j = y 0 + 2a sin(jπ/3).

(

) 2
These are also the translation vectors of length 2a used to generate the tiling. In other words, the centroid coordinates of an innite hexagonal tiling satisfy:

C = {(x 0 , y 0 ) + k(a, √ 3a) + l(2a, 0), ∀k, l ∈ Z}. (3)
As as comparison, the reference square has 8 neighbors that are distributed with an angle of 45°in the square tiling. The edge of square tiles is chosen equal to √ A. With an area of 9A, the Moore neighborhood has more surface than the hexagonal neighborhood. The main dierence is in the distance between neighboring centroids: it varies by a factor √ 2 from even to odd local numbers. This signies that moves in one direction or another may have a dierent weight.

The performance of the hexagonal tiling in the representation of simplied watersheds is evaluated in Section 4.

Tiling and biogeophysical descriptors

Biogeophysical descriptors [START_REF] Grigal | Site and landscape characterization for ecological studies[END_REF] are used in a number of environmental applications. When a hexagonal tiling is applied to partition the DEM and CLC rasters, these scalar descriptors can be assigned to the tiles by averaging and/or by homogenizing the information encoded at the pixel level in the rasters.

Given a DEM and a user-dened outlet, GIS tools such as ArcGIS, QGIS, or TopoToolbox allow to carry out a watershed delineation. The slope and ow direction information are then computed as derivatives by applying some 2D rst order nite dierence operator to the DEM. At the pixel level, such a computation may be implemented by means of the Moore pattern for instance. Then, the resulting ow direction raster is used to deduce the ow accumulation raster and the stream network as a polyline. This process has to be adapted to be suited for a hexagonal tiling.

Two height approximations could be considered for each watershed tile. At a rst glance, the average height seems suitable, however, this choice may create sinks in the tiled watershed model. Thus the minimum height is preferred to conserve the monotonic decrease in elevation from upstream tiles to the outlet tile. The virtual outlet of the tile is xed in the middle of the nearest edge. This becomes the ow face.

A particular attention has to be brought to most upstream tiles and the downstream tile (related to the outlet) because they may partly lie outside the watershed.

For hydrological conceptual models [START_REF] Bergström | Development of a conceptual deterministic rainfall-runo model[END_REF][START_REF] Lindström | Development and testing of the hype (hydrological predictions for the environment) water quality model for dierent spatial scales[END_REF], useful properties at tile level are the areas of the tiles that overlap the watershed, the ow face and the downstream cell, and a set of soil parameters that are partly determined by the land cover.

In this simplied landscape representation, a dominant land cover may be assigned to each tile. In the CLC case, this can be dened using either the CLC nomenclature 1 for a coarse description, or the CLC nomenclature 3 for more features.

Implementation details about biogeophysical descriptors are provided in subsection 3.2.

Tiling and descriptor implementation

An overview of XGeoTiles framework for the tiling of watersheds with hexagons is presented in Fig. 4.

As a preprocessing, elevation and land cover les are loaded into a QGIS project so as to align the DEM and CLC rasters with each other, see Appendix A for details.

The tiling process comprises both the delineation of the watershed from a user dened outlet by means of Topotoolbox, subsection 3.1, and the generation of a hexagonal grid, subsection 3.2. Tiles of the hexagonal grid that overlap or partly overlap with the resulting drainage basin are gathered into the watershed tiling structure denoted by XMap. This process is implemented in the Matlab programming language.

The tiling structure allows for a number of additional prints, subsection 3.3. The 2D boards may be printed as Matlab gures. The hexagonal tiling data information is implemented as an OpenSCAD script to generate a set of joined hexagonal prisms, the height of which accounts for the tile elevation.

DEM analysis in Matlab

The Matlab Mapping toolbox provides very few GIS features when compared to dedicated software. In the Matlab framework, the TopoToolbox library is an object-oriented open source toolbox dedicated to the analysis of DEMs. It provides hydrologists with a much better experience of watershed handling than the Mapping toolbox.

In TopoToolbox, the DEM raster is loaded as an object of the @GRIDobj class, hereafter denoted by DEM. The main properties of the DEM object are for the elevation data stored as a matrix DEM.Z and a 3-by-2 ane transformation as in the Mapping toolbox. This small reference matrix comprises the geographical coordinates of the upper left (reference) pixel of the DEM raster, and the two translation vectors needed to recover the geographical coordinates of the other pixels.

An object of the @FLOWobj class is then built by carving the lled DEM. It imposes a monotonically decreasing sequence of heights along the ow directions and removes sinks. This ow representation allows to deduce the stream network as an object of the @STREAMobj class. A user-dened outlet chosen on the stream allows for the delineation and the storage of the corresponding drainage basin as a GRIDobj object. Such a process underlies the generation of Fig. 1.

As discussed in subsection 2.1, the XGeotiles framework assumes that the land cover information is provided as a raster aligned with the DEM raster. This raster is also provided as a @GRIDobj object. Note that TopoToolbox provides methods for the rasterizaton of shapeles, the reprojection and the resampling of GRIDobj variables. However, to facilitate the uptake by geographs, we recommend to align the two rasters as described in Appendix A.

As shown in Figs. 1 and2, the bounding box of the drainage basin can be used to crop the two rasters for a better rendering and for to manage smaller grid objects in the tiling process.

Tiling structure

The hexagonal partition into tiles of area A is carried out from the reference pixel in the directions indicated in the reference matrix following equation [START_REF] Den Haan | The virtual river game: Gaming using models to collaboratively explore river management complexity[END_REF]. In this formula, bounds for indices k and l may be deduced from the pixel size and the bounding box of the DEM. The watershed tiling is limited to the set of tiles that overlap or partly overlap the drainage basin. A threshold, for instance an overlap larger than 0.5A, may be imposed.

As an illustration, see Fig. 5. The Moderbach's watershed DEM is partitioned into 89 tiles of 1 km 2 , the color of which corresponds to their minimum height. This tiling into hexagon is named XDEM. The stream network computed with the DEM is still plotted to ease the understanding. The construction of a simplied stream network (called XStream) and adapted to the tiling is discussed in subsection 3.2.1. A Matlab structure denoted XMap is chosen to store tile and descriptor information into elds. The naming convention is displayed in Tab. 1 together with their siting (tile centroid, tile or neighborhood (NBH)) and their unit. The tile inradius and circumradius can be deduced from the vertices coordinates. Fields for the mean height or the soil properties may be added to this minimal XMap structure. The XMap structure may be saved as a .csv le using the writetable function of Matlab.

Simplied streamow

For most watersheds, the stream network exhibits a tree-like pattern and sub-basins have a unique outlet. This may not be the case when a tiling is implemented.

The choice is made to assign a unique outlet to the tile, located at the point of the minimum elevation of the DEM pixels inside the tile and the watershed. Thereby, the carving preprocess ensures that this pixel lies near to the tile boundary and determines the outlet face and the downstream tile.

As implemented, this coarse but functional stream network XStream suces to reestablish a monotonic height decreasing sequence of tiles running from the upstream tiles to the outlet tile. The simplied XStream network is plotted by linking one tile after another the centroid to the centroid of its downstream tile. Clearly there are a number of deviations, the geographic extent depends on the tile area and the distance between the hills in the landscape. A rened tiling improves the accuracy of the representation as shown in section 4.

Land cover indicators

In a nutshell, the CLC information is encoded using a three-level hierachical nomenclature that distinguishes 5 main categories at level 1 (nomeclature 1), 15 categories at level 2, and 44 categories at level 3 (nomenclature 3). The CLC codes are strings of 3 digits. The rst digit encodes the main category, see Fig. 2. This allows to simplify the representation and/or to initialize a game board.

For now, the dominant land cover is assigned to the tile. Figure 7 draws the resulting hexagonal tilings XCLC1 and XCLC3, computed with the nomenclature 1 (left picture) and the nomenclature 3 (right picture), respectively. The dierences are mainly due to the mosaic of distinct agricultural practices taking place on a tile. Although a tile may be mostly dedicated to agriculture, the dominant cover may be the urban fabrics as for the central red tile of XCLX3. A similar explanation arises for the blue tile of XCLC3. As can be seen in Fig. 2, the lower urban fabrics represented in XCLC1 divides in two categories in the nomenclature 3, the areas of which are smaller than the area of the waterbody. Again, some renement improves the representation as demonstrated in section 4.

Board rendering and printing

Apart from the image le format proposed in Matlab, the XMap data contain the information needed for the implementation of a number of less conventional plots and prints. Among them, shapeles for QGIS and STL les for 3D printers are the most promising in the context of serious games.

On the one hand, the XMap elds required to generate a land cover shapele are the ID, the vertex co-ordinate arrays X and Y, the CLC codes MainCLC1 and MainCLC3 for nomenclatures 1 and 3 respectively, and an additional Geometry eld set to Polygon for all the tiles. These are gathered in a new Matlab structure to be saved with the shapewrite method. The resulting shapele may be viewed using QGIS and a .qml le adapted to the CLC data. This geographical hexagonal tiling may be used within QGIS to clip any other raster data by extent so as to project the resulting image on a board.

On the other hand, the centroid coordinates X0, Y0 and minimum height Z0 stored in the XMap structure may be used to pave the watershed with vertical hexagonal prisms. Asuming scales of 1:100000 in the horizontal plane and 1:5000 for the heights, the resulting prismatic model may be printed on a small 3D printer. In practice, the list of tiles of XMap, here dened by (X0,Y0,Z0), are turned into a list of calls of the Xprism module reported Appendix B. This is interpreted by OpenSCAD, Fig. 8, to produce a stl le. This small scale prototype of a watershed board involving 89 tiles has been successfully printed on an Original Prusa MINI+. The TikZ package, used for Fig. 3, may also be used to create such graphic hexagonal or prismatic elements in L A T E X.

Experiments and discussion

The reliability of the hexagonal tiling is evaluated with respect to hydrology and land cover representation on rened grids to exhibit strengths, weaknesses and opportunities.

Watershed tiling and hydrology

Conceptual hydrological models such as [START_REF] Bergström | Development of a conceptual deterministic rainfall-runo model[END_REF][START_REF] Lindström | Development and testing of the hype (hydrological predictions for the environment) water quality model for dierent spatial scales[END_REF] rely on a very few soil and vegetation parameters to link the rainfall over a watershed to the discharge at its outlet. Therein, the area acts as a scaling factor applied to calculations operated on a unit of eld. This subsection analyses how the hexagonal tiling complies with a conceptual model by evaluating discrepancies on the total area and by comparing the visual agreement between the DEM and XDEM stream networks. Computations are carried out by using tiles of 0.1 km 2 and a threshold of 0.5 to determine whether a hexagon is or not accounted in the watershed tiling. Table 2 shows some topographical characteristics of tiled watersheds computed for a range of areas. On the one hand, tiles of 1 km 2 may be useful to play a serious game with a relatively small number of tiles. On the other hand, smaller tiles yields more accurate representations for the height. The relative error on area is of the order of 1% in all the three experiments.

Figure 10 plots the DEM stream network and the stream network built during a tiling design process into 10 times more hexagons than used in Fig. 6. As already suggested in subsection 3.2.1, the use of smaller tiles improves the agreement between these two networks.

Land cover tiling and serious game

The water-food-energy security nexus is central to sustainable development. Agriculture is a key actor in all three, consuming freshwater and energy to produce food. Consequently, Utilized Agricultural Land (UAL) may be a key criterion in serious games about water management. Within the proposed tiling, UAL corresponds to tiles with yellow shades. Figure 11 presents the land cover tiling into 897 tiles of 0.1 km 2 using the nomenclatures 1 and 3. Such a renement improves the visual agreement between the land cover tilings and the original land cover map considerably. As expected, smaller hexagonal tiles promote the land cover parcels with a aspect ratio near to 1. This is especially the case for agriculture and forest surfaces. Villages with a compact shape are better represented. The complex shape of the small waterbodies is captured at the raster pixel level only.

Errors between the CLC1 and XCLC1 pixel values are plotted in Fig. 12 with the CLC1 colormap to assess this statement for the coarse (89) and rened (897) tilngs. The errors evaluated as: disc = abs(CLC1.Z -XCLC1.Z); DIF1 = (disc./disc). * CLC1.Z;

(4) take place at the interfaces between land categories. Note that the non arithmetic values (NaN) that result As expected, the error calculation, Fig. 13, carried out for the nomenclature 3 reveals more dierences in agricultural and forested areas because more land categories are considered. Land cover ratios are measured in Tabs. 3 and 4, respectively. For most of the indicators, the smaller the tile area is, the smaller are the deviations to the watershed ratios. In particular, such a trend may be observed for the large indicators (UAL, forests and articial urban area). The very small land cover indicators (water bodies) are more sensitive. Note that, prior to their usage in a serious game, the land cover indicators may be modied in the XMap structure in a simple manner so as to propose a board that could better satisfy the landscape representation by the players. Furthermore, the game may be initialized with the XCLC1 tiling to propose a simplied representation. Then, the XCLC3 tiling may be used to evaluate if and how the players' actions on UAL tiles have changed their landscape management during the game.

Conclusion

The XGeoTile framework is proposed for the representation of actual watersheds as 2D and 3D boards as to preserve the hydrological and geographical meanings. Hexagonal tiles and prisms are favored as they oer more directions for the moves of water and other resources from one tile to the others.

Although serious games on water governance/management are the main intended application, the numerical experiments prove the reliability of the representation with respect to the area, which acts as a scaling factor in any hydrological simulation. Similarly, the land cover information is well represented as soon as a sucient granularity of tiles is considered.

Combined with hydrological simulations, these boards will be tested during workshops involving farmers and stakeholders in a serious game [START_REF] Barreteau | ExpEau: un jeu pour explorer les dynamiques d'interdépendance dans les socio-agro-hydrosystèmes[END_REF] dedicated to the management of the water quality in rural watersheds.

A. Preprocessing data

For the sake of reproducibility, this appendix describes the data collection and preprocessing steps for the Moderbach watershed (Région Grand Est, France).

Data & providers BD ALTI ®, free 75 m raster le, Regional Corine Land Cover, free vector le,

QGIS project

Add Raster Layer: BDALTIV2_75M_FXX_0975_6900_MNT_LAMB93_IGN69, Add Vector Layer: CLC12_RACAL_RGF and dene the SCR of the project using this CLC le. Note the qml le is missing.

Rasterize CLC12_RACAL_RGF as follows:

* Open Attribute • icode12 as eld,

• Georeferenced units, • Width=Height=75,

• Output extent of the DEM, Clip the two rasters by (the same) extent and save them in GeoTIFF format. This CLC rasterization is described with more detail in websie.cefe.cnrs.fr/formations/qgis-2020.

B. OpenSCAD model

The OpenSCAD model is based on translations in the horizontal plane of hexagonal prisms, the circumradius of which is equal to 1. The height of the prisms is variable.

The 3D board is built as the union of prisms generated by the successive calls to the hexa2scad function that creates the map2scad.txt le. This contains calls to the Xprism module. The small code written in listing B allows to generate Fig. 8 and the corresponding stl le.

Figure 1 :

 1 Figure 1: DEM of the Moderbach's catchment delineated using TopoToolbox [9] (extracted from BD Alti ®75m -V2).

Figure 2 :

 2 Figure 2: Land cover of the Moderbach's catchment. Red shades: articial surfaces, yellow shades: agriculture, green shades: forests, blue shades: wetlands and waterbodies (extracted from Corine Land Cover 12.)

Figure 3 :

 3 Figure 3: Tiling patterns and local numbering for the hexagonal and the Moore neighborhoods.

Figure 4 :

 4 Figure 4: 2D and 3D board design process.)

Figure 5 :

 5 Figure 5: Watershed tiling into 89 tiles of 1 km 2 , threshold of 0.5.

Figure 6 :

 6 Figure 6: Comparison of the stream networks for a watershed tiling into 89 hexagons of 1 km 2 .

Figure 6 compares

 6 Figure 6 compares XStream to the DEM stream net-

Figure 7 :

 7 Figure 7: Land cover, 89 tiles of 1 km 2 , nomenclatures 1 and 3.

Figure 8 :

 8 Figure 8: OpenSCAD model for the watershed tiling.

Figure 9 :

 9 Figure 9: Watershed XDEM, 897 tiles of 0.1 km 2 .

Figure 10 :

 10 Figure 10: Comparison of the stream networks for the watershed tiling into 897 hexagons of 0.1 km 2 .

Figure 11 :

 11 Figure 11: Land cover tiling into 897 tiles of .1 km 2 , nomenclatures 1 and 3.

Figure 12 :

 12 Figure 12: Dierences on land cover representation for the partitions into 89 and 897 tiles, nomenclature 1.

Figure 13 :

 13 Figure 13: Errors on land cover representation for the partitions into 89 and 897 tiles, nomenclature 1.

1 :

 1 x ,y ,z ) 9 // adapted from thingiverse / thing :277727 10 // hexagonal prism ; circumradius =1 11 translate ([x ,y ,z /2General OpenSCAD le for Xboard printing. CRediT authorship contribution statement I Charpentier: Conceptualization of this study, Methodology, Software, Writing.

Table 1

 1 Fields of a watershed tiling structure XMap.

	Fields	Siting Units Comment
	ID	Centroid		Unique identier ()
	X0	Centroid m	Abscissa (m)
	Y0	Centroid m	Ordinate (m)
	X	Tile	m	Vertex abcissas
	Y	Tile	m	Vertex ordinates
	Area	Tile	m 2 Of the overlap with the
				watershed (m 2 )
	Z0	Tile	m	Outlet elevation
	OutFace Tile		Outlet face (local numbering)
	DownTile NBH		Downstream tile ID
	MainCLC1 Tile		Main code, nomen. 1
	MainCLC3 Tile		Main code, nomen. 3

Table 2

 2 Watershed tilings at dierent resolutions.

		Tile area (km 2 )	Watershed
	Tiling	1.0	.25	.10
	Threshold	.5	.50	.50
	#Tiles	89	362	897	1
	Area (km 2 )	89.0	90.5	89.7	89.4
	max. height (m) 265.1 278.5 283.3	304.8
	min. height (m) 210.4 209.7 209.8	209.7
	Time (s)	10.3	22.2	45.8	2.5

Table 3

 3 Land cover ratios, nomenclature 1.

			Tile area (km 2 )	Watershed
	Nomen. 1	1.0	.25	.10
	Artif. surf.	.0337 .0608 .0792 .0958
	Agri. area (UAL) .7303 .6851 .6533 .6597
	Forests	.2360 .2376 .2375 .2226
	Wetlands	.0000 .0000 .0021 .0043
	Wat. bodies	.0000 .0165 .0279 .0176

Table 4

 4 Land cover ratios for nomenclature 3, summed with respect to nomenclature 1.

			Tile area (km 2 )	Watershed
	Nomen. 3	1.0	.25	.10
	Artif. surf.	.0337 .0636 .0803 .0958
	Agri. area (UAL) .6854 .6657 .6533 .6597
	Forests	.2697 .2486 .2353 .2226
	Wetlands	.0000 .0028 .0032 .0043
	Wat. bodies	.0112 .0193 .0279 .0176

Table , *

 , Run the eld calculator to convert the CLC CODE_12 string values to integer values,• Create a new eld icode12, with type "integer", • Use Expression: to_int("CODE_12"),

* Rasterize using:

Acknowledgements

The author warmly thanks the researchers of the LTSER-FR Réseau des Zones Ateliers that participate to the transversal action Exp'Eau [1] for precious discussions on serious games and their commitment in planed workshops.


Jean-Baptiste Clément

Damien Sous
email: damien.sous@mio.osupytheas.fr

Frédéric Golay

Mehmet Ersoy

Wave-driven Groundwater Flows in Sandy Beaches: A Richards Equation-based Model

Keywords: Groundwater, swash zone, seepage, forced coupling, numerical modelling

A groundwater model is developed to simulate flows under the swash zone of sandy beach. Variably saturated beach is described by Richards equation, which is solved by the Rivage code thanks to a discontinuous Galerkin method. The SWASH code is used to simulate wave propagation in nearshore waters and compute suitable boundary conditions at the beachface to force Richards equation. An idealized case is considered to investigate the beach groundwater response to the action of swash cycles associated with long IG-like waves: infiltration/ exfiltration, water table, hydraulic head and pore velocities are outlined. A good qualitative agreement is found with observations from experimental studies in the literature.

INTRODUCTION

Groundwater dynamics in sandy beaches is a matter of interest since beach represents interface between land and open sea and accounts for a wide range of processes: biogeochemical cycles, sediment transport, diffusion of dissolved materials, fresh/salt water exchanges between ocean and coastal aquifer, etc. It is therefore directly involved in the anthropogenic issues in nearshore areas in the global context of climate change and growing urbanization: beach erosion, submersion, aquifer salinization, etc. Many laboratory and field studies have been devoted to understand groundwater flows inside sedimentary beaches (Heiss et al., 2015;Steenhauer et al., 2011). Observational approaches are however still limited to resolve the full spatio-temporal variability of governing variables, in particular the fields of head, saturation and flux. Numerical modelling stands out as a promising prospect, allowing to evaluate more accurately water table fluctuations, flow velocities or infiltrations/exfiltrations during swash events which are for now very difficult to measure directly. Nevertheless, up to now, few models have been devoted to the simulation of wavedriven groundwater flows and their scope remains limited (Li and Barry, 2000;Malott, O'Carroll, and Robinson, 2016).

The swash zone is a dynamic region subject to wave action. Waves are covering and uncovering alternately the beach whose boundary is mobile, porous and characterized by complex flows of various natures which affect groundwater. That is another reason why grasping its dynamics as a whole is difficult. Laboratory experiments and field studies have revealed circulation pattern during the swash cycle. The flow is mainly directed offshore, with inflow in the upper part of the swash zone and exfiltration in the lower part (Sous et al., 2013;Sous et al., 2016;Turner et al., 2016). However, such features are expected to strongly depend on the saturation state of the beach and the duration and location of uprush/backwash phase, with significant consequences on seepage flows or rise of water table (Turner and Masselink, 1998;Turner and Nielsen, 1997).

This paper introduces a groundwater model of Rivage code, developed to provide insight into the circulation under sandy beach swash zone. A key-point of the study is the choice of Richards equation usually used to simulate water flow within unsaturated porous media. Due to the various space and time scales involved by the simulation as well as the strong non-linearities of Richards equation, the numerical scheme needs to be robust and versatile (Farthing and Ogden, 2017). The solution is performed through a discontinuous Galerkin framework alongside time adaptation. This method is currently a high trend in numerical simulation to reach desired resolution and accuracy (Dolejší and Feistauer, 2015).

In the present study, an idealized long wave swash event is tested to evaluate the ability of the model to simulate groundwater response to swash cycles. To this end, a weak coupling is considered: wave-resolving bed pressure data provided by the open-source SWASH code (Zijlema, Stelling, and Smit, 2011) is used as a beachface boundary condition in the Rivage code. Qualitative comparisons are carried out between numerical results from the model and observations from laboratory and experimental measurements found in the literature. Distinctive features of the unconfined coastal aquifer during swash cycles are selected and commented. In particular, discussion focuses evolution of seepage location across the beachface as well as global and local pattern of Journal of Coastal Research, Special Issue No. 95, 2020 pressure and velocity field. The goal of this paper is to present a first case handled by a promising numerical tool. The findings will assist future studies to quantify effects of swash on groundwater flow in sandy beaches for more realistic cases.

METHODS

A numerical model is developed to simulate wave-driven groundwater flows in sandy beaches. For this purpose, Richards equation is used together with boundary conditions monitored by external data. Emphasis is given to Richards equation model.

Groundwater Flow Model: Richards Equation

Richards equation is a classic nonlinear parabolic equation to describe flow in variably saturated porous media:

( ) ( ) ( ) ( ) ( ) 0 t K z θ ψ ψ ψ ∂ -∇ ⋅ ∇ + = (1)
where, ψ denotes the pressure head (cm), K the tensor of hydraulic conductivity (cm/h), z the elevation (cm) and θ the water content. Throughout this paper, hydraulic head (cm) is also used: h z ψ = + . Equation ( 1) is called the mixed form of Richards equation because it involves both variables, θ and ψ . This form was chosen because it is mass conservative and is applicable to a wide range of porous media (Farthing and Ogden, 2017). To close Equation (1), two constitutive laws are necessary for water content and hydraulic conductivity. Several relations are available. In this paper, analytical expressions given from an experiment (Vauclin et al., 1979) are used to model fine sand of fairly regular grain-size distribution, less than 1 mm diameter:

( ) ( ) ( ) ( ) ; , 0 ; ; , 0 s I K K s s B D s K A C if A C if θ ψ ψ θ ψ ψ ψ θ ψ     ≤    =  + +    >   (2)
where, A , B . , C , D are empirical parameters, s K is hydraulic conductivity at saturation and s θ is water content at saturation. It can be noticed that unsaturated and saturated zones are defined according to negative and positive ψ values respectively. Reference (atmospheric) pressure is set as 0 and locates water table. Negative values of ψ can be seen as suction where capillary forces are acting. Effective saturation e S is the water content which has been normalized. Flux (cm/h) is computed as:

( ) ( ) z ψ ψ = - ∇ + q K (3)
Richards equation is a widely used model for water flow in unsaturated soils and rocks (Szymkiewicz, 2013) and thus will be a priori suitable to model groundwater flow in sandy beaches. However, its limitations have to be kept in mind. First of all, flux is modelled as a generalized Darcy's law for two-phase flow, sometimes referred as Darcy-Buckingham's law. This is an empiric law which was first derived for creeping one-phase flow. Two limitations can be expected. First, this model is valid only where advective inertial forces are small compared with viscous forces. This is a strong assumption to be checked since turbulent and fast dynamics can occur at the beach boundary. Nevertheless, this study considers sandy beaches which act like low-pass filter. Small enough velocities are rapidly found beneath the boundary (Sous et al., 2016). The second point is generalization from onephase to two-phase flow which lacks rigour. Finally, the wellknown main hypothesis is elimination of air equation from the diphasic system to keep only water equation, which becomes Richards equation. Indeed, air viscosity is about 55 times smaller than the water viscosity. Consequently, it is assumed that pressure gradients balance faster in air than in water phase, and if air phase is connected continuously at every point with the atmosphere, then air equation is neglected. It means features like trapped air pockets or phenomena driven by air phase cannot be modelled (Szymkiewicz, 2013).

In regards to swash zone, Richards equation is a good, but still poorly used, candidate to explore both saturated and unsaturated beach dynamics. In particular, extensive experience is available in the subsurface hydrology community especially for its numerical resolution and its coupling with surface flows.

Richards equation holds strong non-linearities due to constitutive laws and boundary conditions. Moreover, it degenerates into an elliptic equation under complete saturation. Sharp wetting fronts can also happen according to initial and boundary conditions. For these reasons, designing a robust and general numerical tool can be troublesome (Farthing and Ogden, 2017). Rivage code was implemented to tackle these challenges. Richards equation is solved by fixed-point iteration with an implicit scheme and time adaptation. The discretization is performed with a discontinuous Galerkin framework to capture the range of scales arising in the swash situation thanks to space and order adaptation (Dolejší and Feistauer, 2015).

Boundary Conditions for Forced Coupling and Seepage

Surface and subsurface flows are a unified system by nature. For technical reasons (space/time scales, computational solvability), they are often divided into two submodels which work together though a coupling. Different coupling schemes are possible (Furman, 2008). The simplest one has been chosen for this study. It is called "forced coupling", "external coupling" or "weak coupling". Based on the solution of surface water submodel, a boundary condition is prescribed for the groundwater model. It means submodels are solved in succession, and so, the surface submodel cannot be corrected by feedback. One benefit of forced coupling is that implementation is straightforward because no infiltration/exfiltration model is necessary. This technic can approximate sandy beaches where pore velocities are too slow to affect surface water. Nevertheless, it goes without saying that coupling should gain in maturity by reckoning with an infiltration model or mass conservative principle, without considering retroaction.

The SWASH code simulates propagation of free surface waves according to selected input parameters at the offshore boundary. It is used here in a two-layer configuration, all parameters being kept by default (Zijlema, Stelling, and Smit, 2011). The space and time resolved bed pressure are extracted in order to enforce Richards equation by a space and time dependant Dirichlet boundary condition for hydraulic head. To complete this approach, a so-called seepage boundary condition (Scudeler et al., 2013) is prescribed wherever there are dry conditions, beyond the swash tip. In this way, the model can capture seepage faces, between Journal of Coastal Research, Special Issue No. 95, 2020 surface water and a higher phreatic surface. The seepage boundary condition reads as follow:

0, 0 0 0, if and q n q n otherwise ψ ψ = ≥ ⋅ ≥   ⋅ =  (4)
In other words, if porous medium is saturated and there is outflow, then water pours out at atmospheric pressure. This condition is also called Signorini's type. It can be seen as a nonlinear Robin boundary condition since the length of seepage face is unknown a priori. The seaward seepage flow along beachface, as water is discharging, does not influence surface water.

Description of the Idealized Beach Case

Rivage code was previously validated with several test-cases for Richards equation. They involve various processes such as infiltration into soil columns or recharge/depletion/dewatering of water table. The aim of this paper is to assess an idealized beach case, forced by bed pressure data from the SWASH code.

The beach considered is a straight slope of ratio 1:35 which spans over 105 meters. Figure 1 shows the beach profile and the mesh used. Y-axis is multiplied by 20 to make it more visible.

Hydraulic properties for constitutive laws of Equation ( 2) are written in the Table 1. Initial conditions are hydrostatic with water table level at same height of surface water: 0 1.5 h = m (Figure 1).

The bottom is an impervious boundary:

0 =  q n
. The water table and head are kept fixed on the inland boundary:

1.5 h = m. The beachface is monitored by the SWACH code in the way described in the precedent section.

Table 1. Parameters of constitutive laws for the idealized beach case, (Vauclin et al., 1979). In SWASH, a sinusoidal time-series of free surface elevation is imposed at the offshore boundary to simulate a long infragravitylike wave of period 50 s and amplitude 0.3 m. These types of waves have been observed to be a major driver of beach groundwater dynamics of the Rousty microtidal dissipative beach (Sous et al., 2016). This results in two swash events on the beachface. Time simulated lasts 120 s.

Parameters Hydraulic conductivity

RESULTS

An overview of groundwater features during a swash cycle is presented. Then, results focus on uprush and backwash events. Some selected snapshots from the simulation show interesting observations for the idealized beach case. Surface water and water table are depicted by blue and white solid lines respectively. Hydraulic head contours are displayed together with velocities vectors.

Overall Pattern for Swash Cycles

A first observation concerns wave propagation. When a wave is moving, groundwater flows in the saturated zone respond immediately. Hydraulic head spreads vertically quite sharply under the wave front. Infiltration/exfiltration occurs at the back/ front of the bore respectively. It is particularly visible in the Figure 2 even if this phase is completely ideal because it's about the first wave moving through hydrostatic conditions.

Throughout the simulation, a nearly hydrostatic state is observed in the upper beach area. Moreover, flow divides into seaward and landward components near the uprush limit. Velocities magnitude is rather small everywhere, about 1 cm/h, except for uprush infiltration at the top where it reaches up to 160 cm/h. Around exit point of water table, a local overpressure is observed. It is probably caused by numerical artefacts due to the lack of mesh resolution for the seepage boundary condition.

An overall pattern can be depicted for swash cycles. Two recirculation zones are attached respectively to the back and the front of the swash tongue. These two groundwater cells are interacting together and are evolving during uprush and backwash, growing or diminishing. Sometimes one vanishes completely; sometimes another one gets trapped by recharge induced by water table elevation.

Pattern during Uprush Event

During uprush, there is downward infiltration while the swash tongue is covering the beach. It is mainly concentrated in the upper zone, through the boundary in contact with the fully saturated part of the beach. Infiltration induces a rise of the water table near the exit point, with, as seen in Figure 3.

Pattern during Backwash Event

During backwash, a low hydraulic head zone emerges at the base of the swash zone where exfiltration occurs. Velocities are mainly seaward under the swash tongue moving back. Meanwhile, a small region of high hydraulic head persists in the corner of the saturated zone where water table have been raised by the precedent uprush. It results in recharge and, as water is percolating inside the porous medium, flow divides between seaward and landward direction. Moreover, above the unsaturated zone which has been covered by uprush, water is infiltrating downward in a sharp wetting front. While the swash toe retreating to the sea, a seepage face is appearing drops, locally but quickly, the water table what gives form to a small hump-like structure? Figure 4 presents these observations.

Pattern between Backwash and Uprush Events

The next incoming wave arrives just after the backwash. It pushes the region of low hydraulic head up. As this zone swept away, it reaches the remaining small region of high hydraulic head induced by the water table elevation and makes it disappear temporally.

The incoming wave is bringing a new high hydraulic head front. This causes infiltration which arises water table again, ending seepage and erasing the hump-like structure of water table. Then, a new cycle begins. This intermediate state is shown in Figure 5.

DISCUSSION

In the simulation, the offshore water level is similar to the inland water table, therefore no steady pressure gradient is observed across the beachface by contrast to field or laboratory studies where cross-beach gradients can be observed depending on boundary conditions (Sous et al., 2016;Turner et al., 2016).

Within the beach, overall velocity magnitude is about 1 cm/h which is fairly comparable to measurements made in Sous et al. (2016). The only exception is for infiltration/exfiltration where maximum velocity order is about 100 cm/h. This seems to be less compared with the experiments of Turner et al. (2016). Further tests are now engaged in the exact configurations of the laboratory (Turner et al., 2016) and field (Sous et al., 2016) experiments to allow direct comparisons.

As described by Perera et al. (2019), the moving bore is carrying a sharp hydraulic front. Exfiltration occurs at its front and infiltration at its back when it meets the groundwater conditions of the precedent swash.

A local superelevation of the water table has also been revealed in these observational studies due to the water table charge by short or IG wave forcing under the swash zone. Perera et al. (2019) showed a domed water table in their study. The present benchmark case will be extended to a longer duration to test if the model is able to represent such feature.

Water table elevation induces charge in this zone as long as the low-pressure zone does not affect it with the ongoing uprush. As a consequence, infiltration with flow division occurs near the uprush limit. Laboratory experiments and numerical simulations (Li and Barry, 2000;Turner et al., 2016) both predict this feature but Sous et al. (2016) do not observe it. During the simulation, seepage face has small length because water is able to pour out from sand quite fast following surface water. As a consequence, a small hump-like structure appears during backwash because of seepage along the boundary. This last feature is not described in any study found, but further tests are ongoing with more realistic configurations and longer wave forcings to better understand involved processes.

The two circulation zones with infiltration/exfiltration in the upper/lower part of the beach are also identified in the simulation from Li and Barry (2000) and experiments from Turner et al. (2016). Infiltration and exfiltration involve greater pore velocities than in other parts of the beach. Like it is depicted by Perera (2019), gradients brought by the wave are much larger than the average values in the beach due to larger differences in water depth across the wave.

CONCLUSIONS

This paper deals with a numerical model supported by the Rivage code and built around Richards equation so that it can simulate wave-driven groundwater flows in sandy beaches. By generating waves, the SWASH code computes pressure on the bed to enforce Richards equation with a Dirichlet boundary condition. For this forced coupling, groundwater does not hold sway over surface water but this influence can be neglected if pore velocity is slow compared to wave velocities as it happens for beaches with fine grains. What is more, exfiltration by seepage can be modelled through the seepage boundary condition prescribed on the beachface wherever there is no enforced pressure.

Based on the simulation results from an idealized beach case, the following outlines can be drawn for a swash cycle:

• Groundwater pressure stays nearly hydrostatic under dry area, in the upper beach. Incoming waves propagate groundwater pressure front.

• The swash tongue causes overall infiltration during the bed flooding what elevates water table next to beachface and creates sharp wetting front in the unsaturated zone lapped by swash. For charge induced by water table, flow divides into landward and seaward components.

• During backwash, seaward groundwater can be observed under the swash tongue moving back to the sea. Water table exit point follows the swash toe with late because of sand hydraulic properties, resulting in seepage.

• Two localized groundwater circulation zones evolve and interact during the simulation. One, associated with low pressure and exfiltration, is attached to the swash base during backwash and to the forefront of the next wave during uprush. The other one, associated with high pressure and infiltration, is attached to water table recharge during backwash and to the swash wave during uprush.

This general scheme is consistent with observations made by experimental studies in the literature. Further investigations are needed to assess, on the one hand, numerical effects on the solution, such as discretization level or boundary conditions and, on the other, physical parameters such as constitutive laws, beach topography or incoming waves.

Figure 1 .

 1 Figure 1. Topography of the idealized beach and mesh used for the simulation.

Figure 2 .

 2 Figure 2. Snapshot of hydraulic head contours at t = 12.5 s during the first wave propagation on calm waters.

Figure 3 .

 3 Figure3. Snapshot of hydraulic head contours at t = 37.5 s for the end of the first uprush. Velocities range from about 0.3 to 0.6 cm/h in the beach. Maximum velocity is 35 cm/h for infiltration.

Figure 4 .

 4 Figure 4. Snapshot of hydraulic head contours at t = 57.5 s for the end of the first backwash. Most velocities are less than 1 cm/h.

Figure 5 .

 5 Figure 5. Snapshot of hydraulic head contours at t = 67.5 s for the next incoming uprush. Most velocities are less than 0.8 cm/h.

ACKNOWLEDGMENTS

This work has been supported financially by Provence-Alpes-Côte d'Azur region (France).

LITERATURE CITED

 Dolejší, V. andFeistauer, M., 2015. Discontinuous Galerkin Method. Springer International Publishing. Farthing, M.W. andOgden, F.L., 2017


Damien Sous
email: damien.sous@mio.osupytheas.fr

Samuel Meulé

Frederic Bouchette

Frédéric Bouchette

What on Earth?!

published or not. The documents may come

A fundamental issue of long-term surveys in sandy nearshore areas is to maintain instruments at fixed positions and, incidentally, recording reliable data. This becomes an enduring challenge when trying to track storms. Recording storms is always a very gratifying job, as soon as instruments have been safely retrieved and the recovered data is processed in a warm, stable and cozy desk. A bit more depressing is the number of lost instruments, buried under meters of sand accumulation, destroyed by rocks, sliced by boat propellers, snatched by stainless frame lovers or just disappeared without a trace... GLADYS, a research group on nearshore hydro-morphodynamics in the South of France (Le Grau du Roi, www.gladys-littoral.org) has long been involved in such nearshore in-situ deployments. After years spent deploying supporting frames weighted by hundreds of kilos of concrete or lead, our strategy is to bury long (3 to 6 m) tubes or thick angle irons. The basic procedure is that the metal profile is first buried using a motorpump to fluidize the sand bed, either underwater or on the dry beach (Figure 1a), such that a sufficient portion remains outside the sand (typically 1 to 2 m) and then the instruments are bolted onto it. This type of structure prevents the occurrence of most issues mentioned above. This nice outcome is true in most cases. A rare exception is when the sand accumulation is so great that the whole metal profile totally disappears. This is precisely what occurred on Rousty Beach during a 2014-2015 three-month winter survey. Rousty Beach is located in the Beauduc Gulf in the national park of Camargue, about 5 km east of the Saintes Maries de la Mer, France. The studied area (43.456589° N, 4.493034 °W) is a nearly alongshore uniform sandy beach. The beach is micro-tidal (tidal range about 0.3 m) with a typical two-bar and gently sloping profile (average slope about 0.01).

Wind can be strong in the Rousty area, frequently reaching more than 15 m/s, both from the southeast and north to north-west sectors. Wind events driving significant fluctuations of the mean water level can be observed, ranging from a setup of 1 m during south-easterly winds and low atmospheric pressure storm conditions to a setdown about 0.5 m during strong offshore winds and associated high pressure conditions. The wave climate is mainly bimodal, with moderate (0.5-1 m) and short (peak periods lower than 6 s) south-western swells and more energetic storm waves (larger than 2 m for 40% of the events) from the south-east associated with onshore winds. A cross-shore network of 12 tubes were deployed from the shoaling zone up to the dune (Figure 1b), to record wave transformation across the foreshore together with groundwater dynamics inside the beach (Petitjean et al., 2016;Sous et al., 2016Sous et al., , 2017)).

After a sudden seaward shift of the berm, the 4 m long tube deployed in the swash zone totally disappeared under the sand. Its position was precisely known by DGPS-RTK measurements, so we knew that the top of the tube was somewhere under 1.5 m of sand. The first part of the search procedure was to dig with a shovel, and then to use the motorpump to try to fluidize enough sand to detect the top of the tube. It is like creating a little pool of shifting sand.

While Samuel Meulé was permanently injecting pressurised water into the sand, Damien Sous was trying to find the tube top, first with the shovel, then with one foot, then two feet, and finally he just went down little by little in this boiling pool of fluidized sand until he had to use his snorkel to breathe. To get a better picture of the conditions: this was winter, we were in wetsuits from the early morning in water about 12°C. The expected miracle finally occurred -Damien was just able to touch the top of the tube with one of his toes. After mutual and understandably warm congratulations, the team started to think a bit deeper. The main issue was to catch the big hole drilled at the top of the tube with a rope and a shackle. Damien spent nearly half an hour trying to perform this operation only with two feet -a lost cause. Without really considering what he was going to do, mindful of the setting sun and losing energy and patience, he decided to dive head first into the fluidized sand pool, just to be able to use his hand to fix the shackle -a question of seconds for sure.

Upside down, holding his breath, his head nearly one meter and a half under the sand bed, trying to re-find this stupid tube in a totally dark, thick and moving fluid, a series of questions suddenly hit his brain. What will happen if the motorpump stops (as it usually stops without any reason several times a day)? Who filled up the gas the last time? And when? A frightening memory came to his mind: when the pump stops in the Rousty swash zone (which is not the same for all kind of sands), only a few seconds are necessary before the sand quickly becomes compact again, as a nice beach sand to build sand castles. Damien could feel some parts of his body experiencing the partial re-compaction of the sand when Samuel pulled a bit away the hose to avoid hurting him. This is precisely the point where he said to himself, 'What on Earth?!' The slightest mishap and he would stay preserved in the sand for quite a while... But the story has a happy ending (Figure 2). After a few seconds more, the shackle was finally caught and Damien escaped from the sandy hell by twisting and wriggling his body. And, of course, he remained smiling when shaking his head totally glued by the sand: 'Just a nice walk in the park for a guy like me,' he said. However, a stiff drink was necessary a bit later to recover a normal state.

Damien still tells this story from time to time, particularly when he sees students and colleagues a bit too confident when using motorpumps.

Figure 1 .

 1 Figure 1. (a) Illustration of burying of instrument-supporting poles by sand fluidization using motorpump at Rousty Beach, France, in December 2014. (b) Experimental setup deployed at Rousty Beach. A cross-shore network of instrumented poles was deployed across the swash zone using sand fluidization, including light angle irons (3 m long) for short-term cabled sensors and heavier (10 cm diameter, 4 to 6 m long, with white and orange stripes) tubes for long term autonomous sensors. (Photos: D. Sous.)

Figure 2 .

 2 Figure 2. The day after the 'What on Earth' event. Thanks to the efforts deployed the day before, the tube has been lifted up to the sand surface, i.e. nearly 1.5 m above its initial position. The game was almost won. (Photo: D. Sous.)