# Robust Inverse Framework using Knowledge-guided Self-Supervised Learning: An
application to Hydrology

Rahul Ghosh  University of Minnesota - Twin CitiesMinneapolisMinnesotaUSA
<EMAIL_ADDRESS>,  Arvind Renganathan  University of Minnesota - Twin
CitiesMinneapolisMinnesotaUSA<EMAIL_ADDRESS>,  Kshitij Tayal  University of
Minnesota - Twin CitiesMinneapolisMinnesotaUSA<EMAIL_ADDRESS>,  Xiang Li
University of Minnesota - Twin CitiesMinneapolisMinnesotaUSA<EMAIL_ADDRESS>,  Ankush Khandelwal  University of Minnesota - Twin
CitiesMinneapolisMinnesotaUSA<EMAIL_ADDRESS>,  Xiaowei Jia  University of
PittsburghPittsburghPennsylvaniaUSA<EMAIL_ADDRESS>,  Christopher Duffy
Penn State UniversityState CollegePennsylvaniaUSA<EMAIL_ADDRESS>,  John Nieber
University of Minnesota - Twin CitiesMinneapolisMinnesotaUSA<EMAIL_ADDRESS>and  Vipin Kumar  University of Minnesota - Twin CitiesMinneapolisMinnesotaUSA
<EMAIL_ADDRESS>
(2022)

###### Abstract.

Machine Learning is beginning to provide state-of-the-art performance in a
range of environmental applications such as streamflow prediction in a
hydrologic basin. However, building accurate broad-scale models for streamflow
remains challenging in practice due to the variability in the dominant
hydrologic processes, which are best captured by sets of process-related basin
characteristics. Existing basin characteristics suffer from noise and
uncertainty, among many other things, which adversely impact model
performance. To tackle the above challenges, in this paper, we propose a novel
Knowledge-guided Self-Supervised Learning (KGSSL) inverse framework to extract
system characteristics from driver and response data. This first-of-its-kind
framework achieves robust performance even when characteristics are corrupted.
We show that KGSSL achieves state-of-the-art results for streamflow modeling
for CAMELS (Catchment Attributes and MEteorology for Large-sample Studies)
which is a widely used hydrology benchmark dataset. Specifically, KGSSL
outperforms other methods by up to 16 % in reconstructing characteristics.
Furthermore, we show that KGSSL is relatively more robust to distortion than
baseline methods, and outperforms the baseline model by 35% when plugging in
KGSSL inferred characteristics.

Self-supervised Learning, Inverse Modeling, Forward Modeling

††journalyear: 2022††copyright: acmcopyright††conference: Proceedings of the
28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining; August
14–18, 2022; Washington, DC, USA.††booktitle: Proceedings of the 28th ACM
SIGKDD Conference on Knowledge Discovery and Data Mining (KDD ’22), August
14–18, 2022, Washington, DC, USA††price: 15.00††isbn:
978-1-4503-9385-0/22/08††doi: 10.1145/3534678.3539448††ccs: Applied computing
Earth and atmospheric sciences††ccs: Computing methodologies Machine learning

##  1\. Introduction

Machine learning (ML) is increasingly being used to solve challenging tasks in
scientific applications such as hydrology, lake sciences and crop yield
monitoring. Consider the case of hydrology, where streamflow prediction is one
important research problem for understanding hydrology cycles, water supply
management, flood mapping, and other operational decisions such as reservoir
release. For a given entity (basin/catchment, we use either term
interchangeably), the response (streamflow) is governed by drivers
(meteorological data e.g., air temperature, precipitation, wind speed) and
complex physical processes specific to each entity (Newman et al., 2006; Bhatt
et al., 2014). These complex physical processes are best captured by the
inherent characteristics of each entity (e.g., slope, land-cover). For
example, for the same amount of precipitation, two basins will have very
different streamflow response values depending on their land-cover type. Such
problems are often solved using a mechanistic forward model $f$ that predicts
the response $y_{t}$, given drivers $\boldsymbol{x_{t}}$ and entity
characteristic $\boldsymbol{z}$. A wide variety of scientific models can be
considered as a mapping function between drivers $\boldsymbol{x_{t}}$ (e.g.
weather drivers, climate forcings), and response $y_{t}$ (e.g., streamflow in
river basin, global average temperature), governed by entity characteristics.
Figure 1(a) shows the diagrammatic representation of this forward model. More
recently, Machine Learning (ML) models (e.g. LSTMs) are able to provide state
of the art performance for many scientific applications (Kratzert et al.,
2019). The reason is that ML models are able to benefit from training data
from a large cross-section of diverse training data and thus can transfer
knowledge across basins. These ML models also learn the mapping from the
meteorological drivers and time-invariant basin characteristics to streamflow
and thus essentially emulate the forward models (Razavi et al., 2013).

(a)

(b)

Figure 1. (a) Forward model which uses meteorological drivers
$(\boldsymbol{X_{i}^{t}})$ and entity characteristics $(\boldsymbol{z_{i}})$
to predict response $(\boldsymbol{y_{i}^{t}})$ (b) The inverse model which
approximates entity characteristics $(\boldsymbol{z_{i}})$ by inverting the
forward process.

In the scientific application of streamflow modeling, these entity
(basin/catchment) characteristics are only surrogate variables of the true
basin characteristics (Beven, 2020) and thus can lead to several challenges.
First, there often exists high uncertainty in hydrological measurement, which
in turn causes corruption in basin characteristics. Uncertainty can arise due
to temporal change, spatial heterogeneity, sufficiency of the characteristic
itself to explain the rainfall/runoff process, measurement error, missing
data, and correlation among characteristics that collectively contribute to
streamflow. Second, the full set of basin characteristics may not be measured
across all the river basins, resulting in the incompleteness of basin
characteristics. Missing characteristics hinder the building of a global model
that can leverage data across multiple basins and constrains the
transferability of models built from one region to another. Finally, some
basin characteristics may be essential in modeling the rainfall-runoff
response relation but may be completely unknown, not well understood, or not
present in the available set of basin characteristics. Thus, the ability to
infer these time-invariant basin characteristics from the time-varying
meteorological and streamflow data is essential for model prediction and
hydrological process understanding.

This paper presents an inverse modeling methodology that can be used to
identify or reconstruct static characteristics of an environmental system
given its input and output over time. Figure 1(b) shows the diagrammatic
representation of this inverse problem. Inverse problems (Ongie et al., 2020)
appears in many fields of engineering when the goal is to recover “hidden”
characteristics of a system from “observed” data. In recent years, deep
learning techniques have shown remarkable success for solving inverse problems
in various fields such as compressed sensing, medical imaging (Senouf et al.,
2019), and many more (see (Ongie et al., 2020) for a recent overview). In
general, the inverse problem is ill-posed, i.e., one may not be able to
uniquely recover the input field given noisy and incomplete observations
(Ongie et al., 2020; Gilton et al., 2021). Moreover, traditional methods for
solving $z$ are compute intensive, as they require a large number of model
runs (especially if $z$ has a large dimension).

In this paper, we propose to compute $\boldsymbol{z}$ given
$\boldsymbol{X^{t}}$ and $\boldsymbol{y^{t}}$ efficiently. Deep learning
methods traditionally solve inverse problems by minimizing a cost function (Ma
et al., 2020) that consists of a data-fit term, which measures how well the
reconstruction matches the observations and a regularizer. These methods,
largely based on convolution operator tend to work for linear inverse problems
such as image denoising, super-resolution, and compressed sensing, but for
capturing time varying physical processes such as ours (non-linear inverse
problem), the traditional method fails. This motivates us to ask the question:
_Can we develop a robust inverse framework that can accurately recover static
features when they are noisy/missing/uncertain through an inverse modeling
process ?_

To answer this question, we design a novel inverse framework leveraging
knowledge from the hydrological domain in a self-supervised learning framework
to implicitly extract complex correlations embedded in the input data. We call
our methodology knowledge-guided self supervised learning (KGSSL). We show
that KGSSL can be used to reduce the uncertainty in the measurement of basin
characteristics, impute basin characteristics and identify unknown basin
characteristics. KGSSL enables the extraction of time-invariant
characteristics autonomously by distinguishing learning between the same basin
to be similar for different years and dissimilar with other basins. In cases
where certain basin characteristics are known, we further add a pseudo-inverse
loss on top of the learned embeddings to guide the learning using the known
basin characteristics.

In the rest of our paper, we demonstrate the usefulness of KGSSL on real
streamflow data (CAMELS (Newman et al., 2015)). Specifically, our method is
used to identify existing and estimate new characteristics of the physical
system. If the system’s characteristics are missing or known but uncertain,
then our methodology can impute them or reduce their uncertainty. We
additionally show its usefulness when for any new basins, we don’t have any
knowledge of static characteristics $\boldsymbol{z}$. We show empirical
results on the three research questions raised in this paper, namely
measurement uncertainty, imputation, and unknown nature of the basin
characteristics. We also present the proposed framework as a potential
solution for performing source specific modulation in the driver-response
model. Our main contributions are listed below:

  * •

We demonstrate the power of leveraging domain knowledge between entities in
the context of an inverse problem.

  * •

We perform extensive experiments on real streamflow data and show that KGSSL
outperforms baseline by 16 % in reconstructing characteristics.

  * •

We further show that our approach is highly robust to uncertainty in the basin
characteristics. KGSSL outperforms the state-of-the-art forward model by 35%
when plugging KGSSL inferred characteristics.

##  2\. Related Works

Inverse problems (Ongie et al., 2020) always exist together with their forward
problem. The goal of the inverse problem is to recover ”hidden” information
(which we cannot observe directly or is very expensive to observe) from
readily available ”observed” data. Unfortunately, the inverse is often both
intractable and ill-posed, since crucial information is lost in the forward
process. However, the inverse process is required to inform us about physical
parameters of the system (e.g., mass, temperature, physical dimensions, or
structure), sources of influence, reconstruction of the coefficients in the
equations that we cannot observe otherwise. Inverse problems are studied for
many environmental science branches, i.e., hydrogeology (Zhou et al., 2014),
geophysics (Kim et al., 2018), oceanography (Woolway et al., 2021),
meteorology (Pecha et al., 2021), remote sensing (Dao et al., 2021), etc. For
example, an inverse problem arises when we reconstruct Earth’s interior by
modeling the physical propagation of seismic waves (Ritsema et al., 2000).
Similarly, in reservoir engineering (Ganachaud et al., 2000), given various
measurements of geophysical fields, an inverse problem arises to determine the
subsurface properties, such as the permeability field. Most of the recent deep
learning approaches (Ardizzone et al., 2018; Asim et al., 2020) model
forward/inverse mapping within a single network. However, in hydrology, the
initial physical parameters are not known reliably (Kumarasamy et al., 2018)
for the basins/catchment due to temporal and spatial heterogeneity. This leads
to a noisy forward operator, which makes existing inverse approaches
ineffective. This motivates us to design a robust inverse framework impervious
to corrupted basin characteristics.

Due to abundant unlabeled data in computer vision, recently, researchers have
started investigating self-supervised methods (Jing et al., 2020) for model
training. In self-supervised learning, the models are trained using pretext
tasks instead of an actual task. For example, image colorization (Larsson et
al., 2017), image inpainting (Pathak et al., 2016), solving image-jigsaw
(Noroozi et al., 2016), learning by counting (Noroozi et al., 2017),
predicting rotations (Gidaris et al., 2018), etc. For a comprehensive
understanding of self-supervised representation learning, we would like to
redirect the reader to a survey by Jing et al. (Jing et al., 2020). Our work
utilizes a self-supervised loss called InfoNCE loss (Chen et al., 2020) to
counter the uncertainty in the basin characteristics by implicitly extracting
complex correlations embedded in the meteorological drivers. The use of
InfoNCE loss in our work is closely related to that of (Ravula et al., 2021),
which trained the teacher-student network using contrastive loss to recover
the true feature of a corrupted image. However, our problem domain (learning
relationship between time-varying complex physical processes) is fundamentally
different from vision-related inverse problems. In addition, our method
differs in the following aspects. First, we employ the self-supervised
learning method in the time-series domain, whereas most of the applications
are in the vision domain. Second, we design novel pretext tasks in hydrology
using domain knowledge. Finally, we focused extensively on robustness in our
work and showcased our methodology’s use on both supervised and unsupervised
learning.

##  3\. Method

In this work, we study the driver-response relation in dynamical systems.
Specifically, we assume a dataset consisting of $N$ entities (an entity can be
a lake, basin or streams in a river-network). For each entity $i$, the daily
drivers are represented by $\boldsymbol{X_{i}}$ as a multivariate time series
for $T$ timestamp i.e.
$\boldsymbol{X_{i}}=[\boldsymbol{x^{1}_{i}},\boldsymbol{x^{2}_{i}},\dots,\boldsymbol{x^{T}_{i}}]$
where $\boldsymbol{x^{t}_{i}}\in\mathbb{R}^{D_{x}}$ indicates input vector at
time $t\in T$ with $D_{x}$ dimension.
$\boldsymbol{z_{i}}\in\mathbb{R}^{D_{z}}$ denotes the static characteristic
vector of an entity with $D_{z}$ dimensions. The transient response
corresponding to ${(\boldsymbol{X_{i}},\boldsymbol{z_{i}})}$ for an entity is
denoted by $\boldsymbol{Y_{i}}=[y^{1}_{i},y^{2}_{i},\dots,y^{T}_{i}]$.

Our proposed method KGSSL, infers time-invariant entity characteristics
($\boldsymbol{z_{i}}$) given the time-varying driver ($\boldsymbol{X_{i}}$)
and response ($\boldsymbol{Y_{i}}$) data. KGSSL has several components. First,
a Sequence Encoder is used to extract a fixed length representation from the
driver-response time-series. Second, a reconstruction loss
($\mathcal{L}_{Rec}$) that forces the fixed length representation to capture
the information stored in driver-response time-series by penalizing bad
driver-response time-series reconstructions. Third, a Knowledge-guided
Contrastive Loss ($\mathcal{L}_{Cont}$) that implicitly extract complex
correlations embedded in the driver-response time-series and enforces the
physical knowledge that the entity characteristics are time-invariant.
Finally, a PseudoInverse loss ($\mathcal{L}_{Inv}$) that encourages robust
reconstruction of entity characteristics from the fixed length representation
using a feed-forward network. Thus, the final loss function for training KGSSL
is

(1)
$\mathcal{L}=\lambda_{1}\mathcal{L}_{Rec}+\lambda_{2}\mathcal{L}_{Cont}+\lambda_{3}\mathcal{L}_{Inv}$

where $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$ are hyper-parameters to
control the weights of three loss terms. $\mathcal{L}_{Inv}$ is added only
when the entity characteristics are known and available for training. We can
also train KGSSL using only the self-supervised loss functions,
$\mathcal{L}_{Rec}$ and $\mathcal{L}_{Cont}$. In the following subsections, we
discuss each of these components in detail and provide intuition behind such
design choices.

KGSSL generates time invariant and entity specific embeddings from driver-
response time-series data. Specifically, for each entity $i$ in a given set of
N entities, we randomly select two sequences of length $W$. Let $S_{a_{i}}$
and $S_{p_{i}}$ be the two sequences taken from the time-windows
$t_{a_{i}}:t_{a_{i}}+W$ and $t_{p_{i}}:t_{p_{i}}+W$, respectively. This
results in $2N$ sequences and each element in these sequences are formed by
concatenating the drivers-response time-series of the entity
($[\boldsymbol{x^{t}_{i}};y^{t}_{i}]$).

Figure 2. Bidirectional LSTM based Sequence Encoder

###  3.1. Sequence Encoder

We use a sequence encoder to encode the temporal information and the
interaction between the driver and response in these sequences. LSTM is
particularly suited for our task where long range temporal dependencies
between driver and response exist as they are designed to avoid exploding and
vanishing gradient problems. However, LSTMs are designed to run forward in
time and cannot provide explainability on the current time-steps given the
future data. To capture this information we use a Bidirectional LSTM based
sequence encoder $\mathscr{E}$ (Figure 2). Specifically, we build two LSTM
structures:the forward LSTM and the backward LSTM. The two LSTM structures are
the same except that the time-series is reversed for the backward LSTM. Each
LSTM uses the following set of equations to generate the embeddings for a
sequence,

(2)
$\begin{split}\boldsymbol{i_{t}}&=\sigma(\boldsymbol{W_{i}}\left[[\boldsymbol{x^{t}};y^{t}];\boldsymbol{h^{t-1}}\right]+\boldsymbol{b_{i}})\\\
\boldsymbol{f_{t}}&=\sigma(\boldsymbol{W_{f}}\left[[\boldsymbol{x^{t}};y^{t}];\boldsymbol{h^{t-1}}\right]+\boldsymbol{b_{f}})\\\
\boldsymbol{g_{t}}&=\sigma(\boldsymbol{W_{g}}\left[[\boldsymbol{x^{t}};y^{t}];\boldsymbol{h^{t-1}}\right]+\boldsymbol{b_{g}})\\\
\boldsymbol{o_{t}}&=\sigma(\boldsymbol{W_{o}}\left[[\boldsymbol{x^{t}};y^{t}];\boldsymbol{h^{t-1}}\right]+\boldsymbol{b_{o}})\\\
\boldsymbol{c_{t}}&=\boldsymbol{f_{t}}\odot\boldsymbol{c_{t-1}}+\boldsymbol{i}\odot\boldsymbol{g_{t}}\\\
\boldsymbol{h_{t}}&=\boldsymbol{o_{t}}\odot\tanh{(\boldsymbol{c_{t}})}\\\
\end{split}$

Each of the forward and backward LSTM takes in a sequence $S$ as input and
generates corresponding embeddings $\boldsymbol{h_{f}}$ and
$\boldsymbol{h_{b}}$ $(\boldsymbol{h}=\mathscr{E}(S))$. These embeddings are
essentially the final hidden states of each LSTM. The embeddings for the
forward LSTM ($\boldsymbol{h_{f}}$) and backward LSTM ($\boldsymbol{h_{b}}$)
are added to get the final embeddings $\boldsymbol{h}$ as shown in Figure 2.
These embeddings capture the temporal information as well as the driver-
response interaction by modeling the change in streamflow due to the weather
drivers in both forward and backward directions.

###  3.2. Reconstruction Loss

To preserve the key information from driver-response data, we use a standard
LSTM based decoder $\mathscr{D}$ that reconstructs the sequence back from the
embedding ($\hat{S}=\mathscr{D}(\boldsymbol{h})$). The LSTM decoder uses its
own output at the previous time-step as the input for the current time-step
and thus can be regarded as a sequence generator using the embedding
$\boldsymbol{h}$ as a prior. The reconstruction error is computed as the mean-
squared error between the reconstructed and the original sequence, as shown
below,

(3)
$\mathcal{L}_{Rec}=\frac{1}{2N}\sum_{e\in\\{a,p\\}}\sum_{i=1}^{N}MSE(\hat{S}_{e_{i}},S_{e_{i}})$

Here $\mathcal{L}_{Rec}$ acts as a regularizer in representation learning, by
extracting meaningful information from the time-varying input data. However,
since we are interested in extracting the time-invariant information from the
time-series data, solely relying on $\mathcal{L}_{Rec}$ leads to sub-optimal
performance. $\mathcal{L}_{Rec}$ promotes preservation of information about
the time-series in the embeddings which will later be used by the decoder to
reconstruct back the input time-series.

###  3.3. Knowledge-guided Contrastive Loss

Each entity’s response to a given driver is governed by complex physical
processes captured by its inherent physical characteristics that remain
constant through time. Moreover, different entities have different responses
to the same driver due to the differences in their inherent characteristics.
We use this physical knowledge of entities to define a self-supervised
contrastive loss (Chen et al., 2020; Oord et al., 2018). Specifically, the
sequences $S_{a_{i}}$ and $S_{p_{i}}$ of an entity form a positive pair, and
for each positive pair, we treat the other 2(N-1) sequences within a batch as
negative examples. Thus, the contrastive loss forces the embeddings
$\boldsymbol{h_{a_{i}}}$ and $\boldsymbol{h_{p_{i}}}$ resulting from the
sequences $S_{a_{i}}$ and $S_{p_{i}}$ of the same entity to be similar and
different from the embeddings of other basins. For a given positive pair, the
loss is calculated as,

(4)
$\begin{split}l(a_{i},p_{i})=&\frac{\exp{(sim(\boldsymbol{h_{a_{i}}},\boldsymbol{h_{p_{i}}})/\tau)}}{\sum_{e\in\\{a,p\\}}\sum_{j=1}^{N}\exp{(sim(\boldsymbol{h_{a_{i}}},\boldsymbol{h_{e_{j}}})/\tau)}}\\\
+&\frac{\exp{(sim(\boldsymbol{h_{p_{i}}},\boldsymbol{h_{a_{i}}})/\tau)}}{\sum_{e\in\\{a,p\\}}\sum_{j=1}^{N}\exp{(sim(\boldsymbol{h_{p_{i}}},\boldsymbol{h_{e_{j}}})/\tau)}}\end{split}$

where,
$sim(\boldsymbol{h_{a_{i}}},\boldsymbol{h_{p_{i}}})=\frac{\boldsymbol{h_{a_{i}}}^{T}\boldsymbol{h_{p_{i}}}}{\|\boldsymbol{h_{a_{i}}}\|\|\boldsymbol{h_{p_{i}}}\|}$.
Thus, the total contrastive loss for 2N such positive pairs is given as,

(5) $\mathcal{L}_{Cont}=\frac{1}{2N}\sum_{i=1}^{N}l(a_{i},p_{i})$

Both $\mathcal{L}_{Cont}$ and $\mathcal{L}_{Rec}$ do not require any
supervised information and thus can work with a large number of entities for
which we only know the driver-response time-series. Moreover, later in the
results, we show that using only one of these losses leads to sub-optimal
performance and, thus we use a combination of these two losses.

Figure 3. The proposed inverse model generates embeddings for a basin from the
LSTM Encoder (Figure 2)and is trained in a self-supervised manner. Strong
supervision $(\mathcal{L}_{Inv})$ is added when ground-truth characteristics
are available for a limited number of basins

###  3.4. PseudoInverse Loss

Reconstruction Loss and Knowledge-guided Contrastive Loss is used to
extrapolate entity characteristics from the time-varying driver
($\boldsymbol{X_{i}}$) and response ($\boldsymbol{Y_{i}}$). However, if some
entity characteristics are known (albeit noisy (Song et al., 2020)/uncertain),
the above loss functions fail to account for them during the model training.
To improve our inverse framework, we propose using PseudoInverse loss that
utilizes incomplete/uncertain missing entity characteristics as a source of
supervision. Specifically, we add a feed-forward layer $\mathscr{I}$ on
sequence encoder output $h$ to estimate
$\boldsymbol{\hat{z}}=\mathscr{I}(\boldsymbol{h})$ and then we define
regression loss with the available set of entity characteristics ($z$) as
shown in Figure 3 . Pseudoinverse loss is defined as follows:

(6)
$\mathcal{L}_{Inv}=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{z}\sum_{j=1}^{z}(z_{i}^{j}-\hat{z}_{i}^{j})^{2}$

###  3.5. Reconstructing static characteristics given temporal data

Our KGSSL framework can be used to generate entity-specific embeddings as well
as static characteristics. Specifically, given input-drivers
$X_{i}=[\boldsymbol{x^{1}_{i}},\boldsymbol{x^{2}_{i}},\dots,\boldsymbol{x^{T}_{i}}]$
where $\boldsymbol{x^{t}_{i}}\in\mathbb{R}^{D_{x}}$ and output-response
$Y_{i}=[y^{1}_{i},y^{2}_{i},\dots,y^{T}_{i}]$ time-series of length $T$ for an
entity, we break the combined time-series into $T/W$ sequences of length $W$.
Each of these sequences $S_{i}^{j}$ are fed to the encoder $\mathscr{E}$ to
generate an embedding $\boldsymbol{h}_{i}^{j}$, which are further fed into the
inverse regressor $\mathscr{I}$ to predict the static characteristics
$\boldsymbol{\hat{z}}_{i}^{j}$. By taking the element-wise mean of the
embeddings, we get the final embeddings of the entity. Similarly, we get the
final estimate of the static characteristics along with their uncertainties by
taking the element-wise mean and standard deviation of the sequence specific
predictions, as shown,

(7)
$\boldsymbol{h}_{i}=\frac{W}{T}\sum_{j=1}^{T/W}\boldsymbol{h}_{i}^{j}\qquad\boldsymbol{\hat{z}}_{i}=\frac{W}{T}\sum_{j=1}^{T/W}\boldsymbol{\hat{z}}_{i}^{j}\qquad\boldsymbol{unc}_{i}=\sqrt{\frac{W}{T}\sum_{j=1}^{T/W}(\boldsymbol{\hat{z}}_{i}^{j}-\boldsymbol{\hat{z}}_{i})^{2}}$

As more and more years of data are made available for an entity, the
embeddings and the predictions of the static characteristics become more
certain and informative.

Figure 4. Experimental setting followed in the paper for training and testing
of the ML models

##  4\. Experimental Results

###  4.1. Datasets and Implementation details

We evaluate KGSSL using the CAMELS (Catchment Attributes and MEteorology for
Large-sample Studies) dataset, which is extensively used for investigating
hydrology processes, in particular, streamflow prediction (Addor et al.,
2017). CAMELS compiles meteorological forcing data (e.g. precipitation, air
temperature), streamflow observation, calibrated physical model simulation,
and catchment characteristics(see Appendix A.1 for a complete list), all of
which makes it possible to leverage recent developments in machine learning,
in particular deep learning, in the hydrology community to advance continental
hydrology modeling (Kratzert et al., 2019; Feng et al., 2020). In particular,
using the CAMELS dataset, Kratzert et al. (Kratzert et al., 2019) showed that
a global scale LSTM model (that uses known static characteristics as input in
addition to weather drivers) can outperform state-of-the-art physics based
hydrological model that are individually caliberated for each basins.

Following the set up used by Kratzert et al. (Kratzert et al., 2019), our
study uses data for 531 basins from CAMELS for the periods (1989-2009). Of
these, (2001-2008) is used for model building, and the rest is used for
testing. Fig. 4 show the experimental setting followed in this paper. Like
Kratzert et al. our study uses 27 basin characteristics organized by
physically meaningful groups: climatology, soils/geologic conditions and
geomorphology/land-cover. These three groups of characteristics can generally
be assumed to represent physical characteristics that contribute more or less
to the rainfall/runoff process in any given catchment.

Figure 5. Scatter plot of the CAMELS Estimates for static characteristics (x-axis) for the testing basins vs. reconstructed annual characteristics (y-axis). The error bar across the points show the variation of the reconstructed characteristics annually across test years. Method | RMSE | CORR
---|---|---
LSTM | 0.540 | 0.795
KGSSL($L_{Rec+Inv}$) | 0.493 | $\mathbf{0.831}$
KGSSL($L_{Cont+Inv}$) | 0.514 | 0.815
KGSSL($L_{Rec+Cont+Inv}$) | $\mathbf{0.465}$ | 0.824
Table 1. Average root mean square error (RMSE) and correlation (CORR) for
$131$ test basins during testing period.

We create input sequences of length 365 using a stride of half the sequence
length, i.e., 183. This results in 13 windows for the data used for model
training and 19 for the testing period. All LSTMs used in our architecture
have one hidden layer with 64 units. The feed-forward network to reconstruct
characteristics has one hidden layer followed by activation to introduce non-
linear combinations of the embeddings. The hyperparameter $\lambda_{1}$,
$\lambda_{2}$, and $\lambda_{3}$ are set at 1, 1, and 1 respectively. The
value of $\lambda_{1}$, $\lambda_{2}$ and $\lambda_{3}$ are selected to
balance the supervised and unsupervised components of the loss function.
Higher values for $\lambda_{3}$ lead to lower training loss but at the expense
of loss of robustness to noise in the static characteristics(see Appendix A.2
for more details about the hyperparameter search). To reduce the randomness
typically expected with network initialization, we train five models with
different initialization of deep learning model weights. The predictions were
then further combined into an ensemble by averaging prediction from these five
models.

In sections 4.2,4.3,4.4 we evaluate the ability of KGSSL to estimate the
entity characteristics in test basins under various conditions, including when
characteristics in the training set are corrupted or missing. Section 4.4
considers the case where the catchment characteristics are not available
during training. In addition, section 4.6 shows the ability of KGSSL to
improve the forward modeling task.

###  4.2. Estimating the entity characteristics

We train our model using 400 train basins and reconstruct the static
characteristics of the remaining 131 test basins. Table 1 reports average root
mean square error (RMSE) and correlation (CORR) for $131$ test basins during
testing period.The entities have different scale values. Since the RMSE value
is not scale-invariant, we also report a correlation metric, which is scale-
independent. Moreover, the RMSE value measures prediction error, whereas
correlation captures the trend.

Group | Indexes |  Sec 4.2 |  Sec 4.3 |  Sec 4.4
---|---|---|---|---
Original |  90% $\mathcal{N}(0,\sigma_{i})$ |  50% $\mathcal{N}(0,2\sigma_{i})$ |  90% $\mathcal{N}(0,2\sigma_{i})$ | 50% Missing | 90% Missing
Climate | $C1-C9$ | 0.935 | 0.906 | 0.890 | 0.854 | 0.933 | 0.870
Soil-Geology | $S1-S10$ | 0.711 | 0.658 | 0.585 | 0.546 | 0.665 | 0.550
Geomorphology-land cover | $G1-G8$ | 0.841 | 0.812 | 0.783 | 0.751 | 0.825 | 0.783
Mean | 0.824 | 0.786 | 0.745 | 0.709 | 0.802 | 0.725
Table 2. The correlation of reconstructed characteristics for test basins in
reference to true characteristics for different levels of noise and missing
values in train basins.

We make the following high-level observations from our results: a) KGSSL,
which uses both supervised and unsupervised loss functions to infer entity
characteristics, has superior performance ($\mathbf{16\%}$ better RMSE) as
compared to LSTM, which was trained using mean square error loss. b)Each of
the self-supervised losses, i.e., $L_{1}$ and $L_{2}$ individually, leads to
sub-optimal performance; thus, combining these two losses with $L_{3}$ helps
capture the complex physical process accurately. Fig. 5 illustrates the
ability of KGSSL to reconstruct the 27 basin characteristics in the CAMELS
dataset. Note that the reconstructed values are annual averages for each year
in the testing period 1989-1999, while vertical lines show the uncertainty
(UNC) in the prediction as described in Eq. 7). Each individual scatter plot
showcases RMSE, Correlation (CORR) and uncertainty (UNC). Note that in general
KGSSL performs well (correlation¿0.8) in 20 out of 27 cases with correlation ¿
0.9 for 14 of them, and of the remaining 7 cases only one has a correlation
¡0.5.

The results (Table 1 and Figure 5) exhibit that KGSSL is able to predict the
characteristics with acceptable accuracy (corr¿0.8) for most of the
characteristics. In general, the average RMSE and average correlation of the
predicted values are 0.465 and 0.824, respectively. However, the
characteristics reconstruction performance varies among the individual
characteristics. The ones with higher reconstruction RMSEs are usually
accompanied with higher standard deviation. The features with satisfactory
reconstruction performance (lower RMSEs and high correlation) are also more
temporally consistent (lower standard deviation). As discussed in the next
paragraph, inconsistencies in reconstruction performance among the individual
characteristics can be reasoned based on domain knowledge and reflect
uncertainties present in the original CAMELS data set. This interpretation of
the modeling results is arguably an important scientific discovery of our
proposed KGSSL framework.

KGSSL inferred all nine climate characteristics quite accurately. This result
is consistent with the fact that the climate characteristics published in the
CAMELS data set are derived directly from the meteorological forcing data
$\boldsymbol{X_{i}}$. We reason that the elev_mean was also quite accurately
inferred (0.132 RMSE and 0.138 standard deviation) because the mean elevation
is related to climate patterns, and this reasoning also holds true for the
catchment slope characteristic, slope_mean, and vegetation characteristics
(gvf_diff, gvf_max, lai_max, and lai_diff, frac_forest), as these should all
be correlated to meteorological characteristics. The remaining seven
characteristics are uncertain by nature because of involved spatial and
temporal heterogeneities. Some of them also possess uncertainties in the
original data source from which they are derived. Most of these remaining
characteristics are soil-related (e.g., carbonate_rocks_frac,
geol_permeability, soil_depth_pelletier) and are derived as spatial averages
from the catchments. Such derivation overly simplifies catchment spatial
heterogeneity, in particular for large catchments. Therefore, this
simplification might explain the large variance recognized in those
characteristics. Furthermore, as mentioned in (Addor et al., 2017), the
spatial gridded data where those soil characteristics are derived are
uncertain and erroneous in certain geographic regions. It also only
characterizes top layer soils and ignores deep soil information. Consequently,
soil related characteristics are poorly constructed ones. In addition,
area_gages2 is the contributing area where surface runoff is generated, and
this is spatially and temporally highly non-uniform due to the spatial
variability of soil properties, spatial variability of antecedent conditions,
and non-uniformity of incident rainfall. Thus, our reconstruction performance
on area_gages2 is also unsatisfactory.

Figure 6. Comparison of the RMSE of the corrupted values (in blue) generated
by adding Gaussian noise ($\mathcal{N}(0,2\sigma_{i})$) to 50% of basins and
the reconstructed values (in green).

(a)

(b)

(c)

Figure 7. Represent pairwise distance matrices for 531 catchments. Fig (a)
Entry (i,j) is the pairwise distance between characteristic vector of
catchments i and catchment j; Fig (b) Entry (i,j) is the pairwise distance
between embedding vectors generated using KGSSL for catchment i and catchment
j (c) Correlation of each dimension of the learned embeddings with each
physical characteristic.

###  4.3. Robustness to Corruption in available characteristics

As highlighted in Sec.1, we expect uncertainty in the characteristics, , and
furthermore the nature of this uncertainty may be due to temporal and/or
spatial variability, lack of representativeness, measurement error and/or
missing data. The inverse model learns generalizable patterns and hence can
potentially denoise the corrupted characteristics. To emulate this uncertainty
in measurement we randomly corrupt 50% and 90% of the characteristics. Three
experimental setting are thus created. First, to 50% of the characteristics, a
Gaussian noise with 0 mean and 2 standard-deviation is added while the
remaining characteristics are left unchanged. Second, to 90% of the
characteristics a Gaussian noise with 0 mean and 1 standard-deviation is
added. Finally, to the same 90% of the characteristics a Gaussian noise with 2
standard-deviation is added. Those scenarios are created to capture two
perspectives: a small number of characteristics can have a high level of
noise, and a large number of characteristics are corrupted with a relatively
low level noise. We train separate models on the training data using the
corrupted values of these three settings, and the basin characteristics were
predicted using the data from the test years for all the basins and compared
to the original values.

We compare the performance of KGSSL trained using the corrupted data with the
KGSSL trained using original catchment characteristics. Table 2 shows the
performance of various methods in terms of correlation of the predicted
characteristics with the original characteristics. The impact of noise on
characteristics varies among groups, which can be explained by their
dependence on weather data that characteristics are learned. Climate
characteristics are the least sensitive ones because their original
characteristics are derived from weather data. Noisy original characteristics
will not downgrade the reconstruction performance because characteristics can
be learned from weather data anyway. Though not directly related to weather
data, geomorphology-land cover characteristics exhibit geomorphology, land
cover patterns that are implicitly characterized from weather data because of
involved plant growing mechanisms and terrestrial processes. Thus, their
reconstruction performance is much less impacted by noise in the training set.
The worst responses in soil-geology characteristics are likely because they
characterize subsurface processes whose interactions with weather data are
relatively negligible. Such limited usable information in weather data for
soil-geology characteristics constrains the capability of our model to learn.
Figure 6 shows the RMSE computed for corrupted (in blue) and reconstructed
characteristics (in green) with respect to true characteristics averaged
across all the 400 train basins for these two models. We can observe that
KGSSL significantly reduces measurement error in characteristics by an average
RMSE of 1.369.

###  4.4. Robustness to Missing characteristics

Representing physical processes, catchment characteristics often serve as a
unique catchment signature. However, owing to the availability of various data
sources, characteristics that represent one region are likely not available in
another region. It therefore creates a common and important application
scenario where a complete set of catchment characteristics across catchments
are not assured. This limitation is more pronounced for cross-continental
catchments whose characteristics are overlapping rather than exactly matching
with each other. For instance, over half of catchment characteristics (e.g.,
main stream length, bulk density) in CAMELS-CN (a version of CAMELS for China)
(Hao et al., 2021) are not included in the characteristics set of the
catchments in the CAMELS dataset being used in this paper (which only contains
basins from USA) (Addor et al., 2017). The same scenario is also present in
CAMELS version for Great Britain (Coxon et al., 2020), Chili (Alvarez-Garreton
et al., 2018), and Brazil (Chagas et al., 2020). In addition, insufficient
understanding of catchment processes will also lead to a select set of
characteristics that miss the opportunity to capture certain hydrological
processes beyond current hydrological understanding. To address this issue,
the KGSSL can potentially estimate catchment characteristics when they are
missing for some catchments. To emulate such a scenario of missing catchment
characteristics, we use a similar set up as in Sec. 4.3. Instead of adding
Gaussian noise, we treated 50% and 90% of the characteristics to be missing.
We trained separate models for each of these settings on train years and train
basins, where $\mathcal{L}_{Inv}$ was calculated and used for training the
model only when characteristics were available. The catchment characteristics
were predicted using the data from the test years for all the basins and
compared to the original catchment characteristics.

The predicted catchment characteristics using data from the test years are
compared to the original catchment characteristics. For the setting with 50%
missing data, the average RMSE and average correlation of the predicted values
are 0.540 and 0.802 respectively, whereas for 90% missing data, the average
RMSE and average correlation of the predicted values are 0.646 and 0.725
respectively. This result suggests that KGSSL can potentially be used to
impute the missing characteristics. Further, Table 2 shows the robustness of
our method, where we predict the characteristics for the 131 test catchments
using the data from the test years using both the models and compare the
prediction performance to the model trained using the clean data (Section
4.2).

###  4.5. Discovering characteristics in the absence of ground truth(known
characteristics)

Here we investigate the ability of KGSSL to identity time invariant
characteristics that may be missing from available characteristics. We train
the inverse model without using any knowledge of available characteristics
that we can use as a constraint. $\mathcal{L}_{Rec}$ and $\mathcal{L}_{Cont}$
are used for training the model using the data from train years for all 531
basins. Further, using the data from the test years the embeddings for each
basin are computed. To empirically demonstrate the characteristics captured by
the learned embeddings, we calculate the pairwise-euclidean distance between
two basins using their 27d physical characteristics (Figure 7(a)) and compare
them with the distances computed using learned embeddings (Figure 7(b)). We
generate (Figure 7(a)) by reordering the rows in the distance matrix computed
using 27d physical characteristics such that basins with the least distances
between themselves are placed close to each other to form a band-like
structure. The exact order of basins used to generate (Figure 7(a)) is then
further applied to the distance matrix computed using learned embeddings to
generate (Figure 7(b)). From the figure we observe similar patterns in both
the distance matrices which shows that KGSSL generates embeddings that
contains meaningful similarity structure between basins. Further, we calculate
the correlation between the learned embeddings with each of the physical
characteristics for 531 basins. Figure 7(c) provides a measure of the relative
contribution of the 3 groups (C1-9 Climate, S1-10 Soils/ Geology, G1-8
Geomorphology/Landcover) for explaining the rainfall-runoff process. The
vertical axis represents the 27 Static Characteristics, and the horizontal
axis is the embeddings ranked from highest average correlation across
characteristics (left) to lowest average correlation (right). Note that S1 and
G3 have a weak correlation across all embeddings, Collectively the Climate
characteristics show the strongest correlation, followed by
geomorphology/landcover. The soil and geology group represent the weakest
correlation. This might be expected since soil , and geologic properties have
high spatial variability as discussed earlier.

###  4.6. Forward Modeling based evaluation

In previous sections, we observed that KGSSL is able to recover
characteristics under missing/uncertain scenarios. In this section, we take
one step further and plug our retrieved values in state-of-the-art
hydrological models to evaluate the gains achieved in streamflow prediction
performance by these retrieved values compared to the missing/uncertain
values.

LSTMs are extensively used for environmental modeling where both static and
time-series variables are supplied as input (here, static characteristics are
repeated at each time-step). However, the original RNN models were not
designed to exploit static data. Recently, EA-LSTM (Kratzert et al., 2019) has
emerged as one of the state-of-the-art ML-based forward models used in
hydrology that processes the time-series meteorological drivers conditioned on
static characteristics. Henceforth, we compare the streamflow prediction
performance of the EA-LSTM model in two settings: KGSSL inferred features and
original basin characteristics. We report Nash–Sutcliffe model efficiency
coefficient (NSE) score for each forward model run.

Method | Mean NSE
---|---
Baseline(uses known characteristics) | 0.704
Baseline (100% missing) | 0.491
KGSSL (10% missing) | 0.697
KGSSL (50% missing) | 0.700
KGSSL (90% missing) | 0.664
KGSSL (100% missing) | 0.669
Table 3. Model performance for different percentages of missing values.

####  4.6.1. Forward modeling with missing entity characteristics

By design, KGSSL is trained using both supervised and unsupervised loss and
has generalizations capability to infer the missing basin characteristics that
can eventually enhance the streamflow prediction when basin characteristics
are missing/not available. We train all models on all 531 basins during the
train years and test the performance during the test years. Table 3 (first two
rows) report performance of state-of-the-art EA-LSTM (baseline) trained with
all and no static characteristics (Kratzert et al., 2019). The baseline model
where all characteristics are present performs $\mathbf{43\%}$ better (mean
NSE) than the baseline model when some or all basin characters are missing.
This shows the importance of the basin characteristics in modulating the
driver-response network.

To evaluate how much inferred basin features help in the forward model, we
randomly treat 10%, 50%, 90%, and 100% of the characters to be missing. We
impute missing characters using our KGSSL pipeline and run it through the
forward model. Table 3 (last 4 rows) report NSE performance when our model was
used to fill in the static characteristics for different percentages of
missing values. We observe that the forward model trained with reconstructed
characteristics from KGSSL with 10% and 50% missing values perform similar to
the baseline trained with all characteristics. Further with 90% missing
characteristics, the forward model observes only 5% drop when compared with
baseline. In addition, for 100% missing characteristics, we use a total
unsupervised setting in our KGSSL framework, i.e., generate embeddings instead
of characteristics. The KGSSL model with no supervision perform
$\mathbf{35\%}$ better to the baseline with 100% missing characteristics. Even
more impressive is the fact that KGSSL with no supervision (last line) is only
slightly worse than the baseline that uses known characteristics. We attribute
this success to our framework’s knowledge-guided component, which implicitly
extracts complex correlations embedded in the input data.

Method | Mean NSE
---|---
Baseline(actual characteristics) | 0.560
Baseline ($0.5\sigma_{i}$ noise) | 0.474
Baseline ($1\sigma_{i}$ noise) | 0.245
KGSSL ($1$ year) | 0.460
KGSSL ($2$ year) | 0.535
KGSSL ($3$ year) | 0.554
KGSSL ($9$ year) | 0.582
Table 4. Forward model performance with corrupted characteristics and using
KGSSL embeddings. In KGSSL (n-year), the n refers to the number of years of
data utilized to learn the embedding. Figure 8. Basin at year 1992 (Best seen
in color)

####  4.6.2. Forward modeling with corrupted entity characteristics

As shown by Kratzert et al. (Kratzert et al., 2019), uncertainty or corruption
in basin-characteristics can have detrimental effect on the forward modeling.
To demonstrate this, we trained the baseline model on the 400 train basins in
the train years and test the performance on the 131 test basins in the test
years. Table 4 (first row) shows the baseline performance of the forward
model. To model uncertainty in the static characteristics, we add Gaussian
noise ($\mathcal{N}(0,0.5\sigma_{i})$,$\mathcal{N}(0,1\sigma_{i})$) to test
characteristics and measure the performance (Table 4 \- $2^{nd}$, $3^{rd}$
row) of the baseline model. As expected, the forward model is susceptible to
noise in the basin characteristic, and the performance drops significantly
with slight noise. Specifically, the mean NSE drops by 50% with a $1$ standard
deviation noise.

If the basin characteristics are corrupted, we can utilize representation
obtained from the KGSSL trained in self-supervised manner using n-years of
observations (note that this approach does not need any information about
characteristics but it does need a small amount of data to create the
embeddings). Table 4 (second set of rows), showcases the power of this
methodology. As expected, the performance improves as we use more data to
generate embeddings. Note that with only 2 years of data, the EALSTM model
using KGSSL(2 year) outperforms the EALSTM using Corrupted
$\mathcal{N}(0,0.5\sigma_{i})$ characteristics. Moreover, KGSSL(9 year)
generated with 9 years of train data outperform the model with actual
characteristics.

In Figure 8 we plot the actual observed streamflow (black dot) as well as the
predicted streamflow using the various settings of the forward model. We
observe that baseline imputation with corruption performs poorly and is
nowhere close to the actual values. We also note that baseline prediction
(blue line) closely matches our KGSSL predicted values using only 3 years of
data (red line). We attribute this good result to our novel pretext task that
is able to handle time varying physical processes.

##  5\. Discussion and Future Work

In this work, we build a novel inverse framework KGSSL, and demonstrate the
power of leveraging domain knowledge between entities in the context of
streamflow. We performed extensive experiments on the hydrological benchmark
dataset and show that KGSSL outperforms baseline significantly by a margin of
16-35 % under various situations. KGSSL is a first-of-its-kind knowledge-
guided framework that implicitly extracts system characteristics given its
driver and response data.

This paper addresses an important problem in the hydrologic domain, which is
societally relevant. For much of the world, where water is kept in basins for
later redistribution for societal needs, additional hydro climatological
forecasting is required, particularly streamflow. Reasonable estimates are
needed of the time patterns of streamflow about each basin and estimates of
the release schedules necessary to meet societal contracts and environmental
laws. As such, KGSSL can yield societally useful products for resource
management and decision making.

We note that the proposed method is general and can add value in other
applications such as computer vision (self-driving car), where additional
features are used to capture variations in light (Dai et al., 2018), weather
(Volk et al., 2019), and object poses (Alcorn et al., 2019).Note that KGSSL
learns representations without focusing on optimizing the response variable
(i.e., streamflow prediction in our hydrology application). KGSSL can be
further extended by combining both the forward and inverse model in a unified
framework that uses the inverse model to generate a representation and then
uses the learned representation to modulate the forward model. Such an
extension can leverage the recent work from task-aware modulation in machine-
learning (Vuorio et al., 2019; Lin et al., 2021; Zintgraf et al., 2019), and
will be considered in future work.

##  6\. Acknowledgement

This work was funded by the NSF HDR Grant 1934721 and NSF FAI Grant 2147195.
Access to computing facilities was provided by the Minnesota Supercomputing
Institute. John Nieber’s effort on this project was partially supported by the
USDA National Institute of Food and Agriculture, Hatch/Multistate Project MN
12-109.

## References

  * (1)
  * Addor et al. (2017) Nans Addor et al. 2017\.  The CAMELS data set: Catchment attributes and meteorology for large-sample studies.  _Hydrology and Earth System Sciences_ 21, 10 (2017), 5293–5313.  https://doi.org/10.5194/hess-21-5293-2017
  * Alcorn et al. (2019) Michael A Alcorn et al. 2019\.  Strike (with) a pose: Neural networks are easily fooled by strange poses of familiar objects.  _CVPR_ (2019).
  * Alvarez-Garreton et al. (2018) Camila Alvarez-Garreton et al. 2018\.  The CAMELS-CL dataset: Catchment attributes and meteorology for large sample studies-Chile dataset.  _Hydrology and Earth System Sciences_ 22, 11 (2018), 5817–5846.  https://doi.org/10.5194/hess-22-5817-2018
  * Ardizzone et al. (2018) Lynton Ardizzone et al. 2018\.  Analyzing inverse problems with invertible neural networks.  _arXiv preprint arXiv:1808.04730_ (2018).
  * Asim et al. (2020) Muhammad Asim et al. 2020\.  Invertible generative models for inverse problems: mitigating representation error and dataset bias. In _International Conference on Machine Learning_. PMLR, 399–409.
  * Beven (2020) Keith Beven. 2020\.  Deep learning, hydrological processes and the uniqueness of place.  _Hydrological Processes_ 34, 16 (2020), 3608–3613.
  * Bhatt et al. (2014) Gopal Bhatt et al. 2014\.  A tightly coupled GIS and distributed hydrologic modeling framework.  _Environmental modelling & software_ 62 (2014), 70–84.
  * Chagas et al. (2020) Vinicius B.P. Chagas et al. 2020\.  CAMELS-BR: Hydrometeorological time series and landscape attributes for 897 catchments in Brazil.  _Earth System Science Data_ 12, 3 (2020), 2075–2096.  https://doi.org/10.5194/essd-12-2075-2020
  * Chen et al. (2020) Ting Chen et al. 2020\.  A simple framework for contrastive learning of visual representations. In _International conference on machine learning_. PMLR, 1597–1607.
  * Coxon et al. (2020) Gemma Coxon et al. 2020\.  CAMELS-GB: hydrometeorological time series and landscape attributes for 671 catchments in Great Britain.  _Earth System Science Data_ 12, 4 (2020), 2459–2483.  https://doi.org/10.5194/essd-12-2459-2020
  * Dai et al. (2018) Dengxin Dai et al. 2018\.  Dark model adaptation: Semantic image segmentation from daytime to nighttime.  _ITSC_ (2018).
  * Dao et al. (2021) Phuong D Dao et al. 2021\.  Improving hyperspectral image segmentation by applying inverse noise weighting and outlier removal for optimal scale selection.  _ISPRS Journal of Photogrammetry and Remote Sensing_ 171 (2021), 348–366.
  * Feng et al. (2020) Dapeng Feng et al. 2020\.  Enhancing streamflow forecast and extracting insights using long-short term memory networks with data integration at continental scales.  _Water Resources Research_ 56, 9 (2020), e2019WR026793.
  * Ganachaud et al. (2000) Alexandre Ganachaud et al. 2000\.  Improved estimates of global ocean circulation, heat transport and mixing from hydrographic data.  _Nature_ 408, 6811 (2000), 453–457.
  * Gidaris et al. (2018) Spyros Gidaris et al. 2018\.  Unsupervised representation learning by predicting image rotations.  _arXiv preprint arXiv:1803.07728_ (2018).
  * Gilton et al. (2021) Davis Gilton et al. 2021\.  Deep equilibrium architectures for inverse problems in imaging.  _IEEE Transactions on Computational Imaging_ 7 (2021), 1123–1133.
  * Hao et al. (2021) Zhen Hao et al. 2021\.  CCAM: China Catchment Attributes and Meteorology dataset.  _Earth System Science Data_ 13, 12 (2021), 5591–5616.  https://doi.org/10.5194/essd-13-5591-2021
  * Jing et al. (2020) Longlong Jing et al. 2020\.  Self-supervised visual feature learning with deep neural networks: A survey.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ (2020).
  * Kim et al. (2018) Yuji Kim et al. 2018\.  Geophysical inversion versus machine learning in inverse problems.  _The Leading Edge_ 37, 12 (2018), 894–901.
  * Kratzert et al. (2019) Frederik Kratzert et al. 2019\.  Towards learning universal, regional, and local hydrological behaviors via machine learning applied to large-sample datasets.  _Hydrology and Earth System Sciences_ 23, 12 (2019), 5089–5110.
  * Kumarasamy et al. (2018) Karthik Kumarasamy et al. 2018\.  Calibration parameter selection and watershed hydrology model evaluation in time and frequency domains.  _Water_ 10, 6 (2018), 710.
  * Larsson et al. (2017) Gustav Larsson et al. 2017\.  Colorization as a proxy task for visual understanding. In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_. 6874–6883.
  * Lin et al. (2021) Xixun Lin et al. 2021\.  Task-adaptive neural process for user cold-start recommendation. In _Proceedings of the Web Conference 2021_. 1306–1316.
  * Ma et al. (2020) Wei-Chiu Ma et al. 2020\.  Deep feedback inverse problem solver. In _European Conference on Computer Vision_. Springer, 229–246.
  * Newman et al. (2015) Andrew J Newman et al. 2015\.  Gridded ensemble precipitation and temperature estimates for the contiguous United States.  _Journal of Hydrometeorology_ 16, 6 (2015), 2481–2500.
  * Newman et al. (2006) Brent D Newman et al. 2006\.  Ecohydrology of water-limited environments: A scientific vision.  _Water resources research_ 42, 6 (2006).
  * Noroozi et al. (2016) Mehdi Noroozi et al. 2016\.  Unsupervised learning of visual representations by solving jigsaw puzzles. In _European Conference on Computer Vision_. Springer, 69–84.
  * Noroozi et al. (2017) Mehdi Noroozi et al. 2017\.  Representation learning by learning to count. In _Proceedings of the IEEE International Conference on Computer Vision_. 5898–5906.
  * Ongie et al. (2020) Gregory Ongie et al. 2020\.  Deep learning techniques for inverse problems in imaging.  _IEEE Journal on Selected Areas in Information Theory_ 1, 1 (2020), 39–56.
  * Oord et al. (2018) Aaron van den Oord et al. 2018\.  Representation learning with contrastive predictive coding.  _arXiv preprint arXiv:1807.03748_ (2018).
  * Pathak et al. (2016) Deepak Pathak et al. 2016\.  Context encoders: Feature learning by inpainting. In _Proceedings of the IEEE conference on computer vision and pattern recognition_. 2536–2544.
  * Pecha et al. (2021) Petr Pecha et al. 2021\.  Determination of radiological background fields designated for inverse modelling during atypical low wind speed meteorological episode.  _Atmospheric Environment_ 246 (2021), 118105.
  * Ravula et al. (2021) Sriram Ravula et al. 2021\.  Inverse Problems Leveraging Pre-trained Contrastive Representations.  _Advances in Neural Information Processing Systems_ 34 (2021).
  * Razavi et al. (2013) Tara Razavi et al. 2013\.  Streamflow prediction in ungauged basins: review of regionalization methods.  _Journal of hydrologic engineering_ 18, 8 (2013), 958–975.
  * Ritsema et al. (2000) Jeroen Ritsema et al. 2000\.  Seismic imaging of structural heterogeneity in Earth’s mantle: evidence for large-scale mantle flow.  _Science Progress (1933-)_ (2000), 243–259.
  * Senouf et al. (2019) Ortal Senouf et al. 2019\.  Self-supervised learning of inverse problem solvers in medical imaging.  In _Domain adaptation and representation transfer and medical image learning with less labels and imperfect data_. Springer, 111–119.
  * Song et al. (2020) Hwanjun Song et al. 2020\.  Learning from noisy labels with deep neural networks: A survey.  _arXiv preprint arXiv:2007.08199_ (2020).
  * Volk et al. (2019) Georg Volk et al. 2019\.  Towards robust CNN-based object detection through augmentation with synthetic rain variations.  _ITSC_ (2019).
  * Vuorio et al. (2019) Risto Vuorio et al. 2019\.  Multimodal model-agnostic meta-learning via task-aware modulation.  _arXiv preprint arXiv:1910.13616_ (2019).
  * Woolway et al. (2021) R Iestyn Woolway et al. 2021\.  Winter inverse lake stratification under historic and future climate change.  _Limnology and Oceanography Letters_ (2021).
  * Zhou et al. (2014) Haiyan Zhou et al. 2014\.  Inverse methods in hydrogeology: Evolution and recent trends.  _Advances in Water Resources_ 63 (2014), 22–37.
  * Zintgraf et al. (2019) Luisa Zintgraf et al. 2019\.  Fast context adaptation via meta-learning. In _ICML_. PMLR, 7693–7702.

##  Appendix A Appendix

###  A.1. Index of Catchment attributes used in our paper

Group | index | Name
---|---|---
Climate | C1 | p mean
| C2 | pet mean
| C3 | p seasonality
| C4 | frac snow
| C5 | aridity
| C6 | high prec freq
| C7 | high prec dur
| C8 | low prec freq
| C9 | low prec dur
Soil geology | S1 | carbonate rocks frac
| S2 | geol permeability
| S3 | soil depth pelletier
| S4 | soil depth statsgo
| S5 | soil porosity
| S6 | soil conductivity
| S7 | max water content
| S8 | sand frac
| S9 | silt frac
| S10 | clay frac
Geomorphology | G1 | elev mean
| G2 | slope mean
| G3 | area gages2
| G4 | frac forest
| G5 | lai max
| G6 | lai diff
| G7 | gvf max
| G8 | gvf diff
Table 5. Table of catchment characteristics used in this experiment,
description of the characteristics is available and defined in (Addor et al.,
2017)

###  A.2. Hyperparameter Tuning

We used grid search over a range of parameter values to find the
hyperparameters, i.e., $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$,
batch_size, the temperature of contrastive loss, learning rate, dimension of
embedding. Specifically, we considered the following possible parameter
values.

  * •

Dimension of embedding: 32, 64, 128,256

  * •

Learning_rate : 0.0005, 0.001, 0.003, 0.005, 0.05

  * •

$\lambda_{1}$ :0.01, 0.1, 1, 10

  * •

$\lambda_{2}$ :1

  * •

$\lambda_{3}$ : 0.1, 1, 10

  * •

Batch_size: 100, 200

  * •

Temp : 0.1, 0.5, 0.7, 1

We trained our model on train basins in the training period and tested the
model in the testing period on testing basins for each combination in the
parameter set. We chose the parameter set with the least avg RMSE in the train
basins during the validation years as the final parameter configuration

###  A.3. Reproducibility

CAMELS input data is freely available on the homepage of the NCAR at
https://ral.ucar.edu/solutions/products/camels. The code is available at
https://tinyurl.com/bdny7fk6 (KGSSL Code Link).


# Experimental demonstration of negative refraction with 3D locally resonant
acoustic metafluids

Benoit Tallon  Univ. Bordeaux, CNRS, Bordeaux INP, ENSAM, I2M, UMR 5295,
F-33405, Talence, France  Artem Kovalenko  Univ. Bordeaux, CNRS, CRPP,
F-33600, Pessac, France  Olivier Poncelet  Univ. Bordeaux, CNRS, Bordeaux INP,
ENSAM, I2M, UMR 5295, F-33405, Talence, France  Christophe Aristégui  Univ.
Bordeaux, CNRS, Bordeaux INP, ENSAM, I2M, UMR 5295, F-33405, Talence, France
Olivier Mondain-Monval  Univ. Bordeaux, CNRS, CRPP, F-33600, Pessac, France
Thomas Brunet  Univ. Bordeaux, CNRS, Bordeaux INP, ENSAM, I2M, UMR 5295,
F-33405, Talence, France<EMAIL_ADDRESS>
###### Abstract

Negative refraction of acoustic waves is demonstrated through underwater
experiments conducted at ultrasonic frequencies on a 3D locally resonant
acoustic metafluid made of soft porous silicone-rubber micro-beads suspended
in a yield-stress fluid. By measuring the refracted angle of the acoustic beam
transmitted through this metafluid shaped as a prism, we determine the
acoustic index to water according to Snell’s law. These experimental data are
then compared with an excellent agreement to calculations performed in the
framework of Multiple Scattering Theory showing that the emergence of negative
refraction depends on the volume fraction $\Phi$ of the resonant micro-beads.
For diluted metafluid ($\Phi=3\%$), only positive refraction occurs whereas
negative refraction is demonstrated over a broad frequency band with
concentrated metafluid ($\Phi=17\%$).

## Introduction

Since the pioneering works reported by Liu et al.[1], locally resonant
acoustic metamaterials have been attracting great attention [2]. One of the
challenging issues has been the achievement of acoustic metamaterials with a
negative refractive index that offer new possibilities for acoustic imaging
materials and for the control of sound at sub-wavelength scales [3]. In
acoustics, the refractive index $n$ is proportional to $\sqrt{\rho/K}$, where
$K$ and $\rho$ are the bulk modulus and the mass density of the material. Many
works have been devoted to the study of double-negative metamaterials [4, 5,
6, 7, 8, 9, 10, 11, 12] for which the two constitutive parameters $K$ and
$\rho$ are simultaneously negative, leading thus to a negative index[13]. It
is worth noting that such a double-negativity condition is not required for
(real) dissipative metamaterials to get a negative index[14]. When they are
non-negligible, losses may play an important role in the effective acoustic
properties of the metamaterials in such a way that the real part of the
acoustic index can be negative for single-negative metamaterial as
demonstrated in underwater ultrasonic experiments[15], and in air at audible
frequencies[16]. Although the latter work reported on the experimental
observation of negative refraction effects within a 2D acoustic superlens,
negative refraction has never been observed with a 3D acoustic metamaterial up
to now.

In that context, we proposed to use a ”soft” approach, combining various soft-
matter techniques, to achieve soft 3D acoustic metamaterials with negative
index composed of resonant porous micro-beads randomly-dispersed in a yield-
stress fluid[17]. By taking benefit from the strong low-frequency Mie-type
(monopolar and dipolar) resonances of these ’ultra-slow’ particles, single-
band[15] and dual-band[18] negative refractive indices were experimentally
demonstrated. The issue of the experimental observation of negative refraction
was then raised for these water-based metamaterials[19], since the energy
attenuation might be significant in these metafluids due to the intrinsic
absorption in the porous micro-beads, and to the strong resonant scattering by
the particles.

In this paper, we report on negative refraction experiments with metafluids,
composed of soft porous silicone-rubber micro-beads, exhibiting a negative
acoustic index at ultrasonic frequencies[15]. In these experiments, the
metafluid is confined in a prism-shaped box with a small angle
($\theta_{\rm{fluid}}=+2$°) and with a very thin plastic surface in order to
be acoustically transparent. As shown in Fig. 1, this metafluid is directly
deposited on a large ultrasonic immersion transducer operating in water over a
broad ultrasonic frequency range (from 15 kHz to 600 kHz). The large
dimensions of this transmitter ensures the generation of quasi-pure plane
waves propagating in the metafluid along the vertical $z$-axis, from the
bottom to the top. Then, the transmitted acoustic beam refracted by the
interface metafluid/water is scanned in the water tank by using a small
ultrasonic probe. Two samples, composed of micro-beads with similar mean
diameters $d$ of about 750 $\mu$m, are considered in this study with two
different volume fractions $\Phi$, referred to as the diluted metafluid
($\Phi=3\%$) and concentrated metafluid ($\Phi=17\%$). The acoustic properties
of the micro-beads are given in a previous work[20]. One of the particular
features of these porous particles is their very low longitudinal phase
velocity $c_{L}$ (= 120 m.s-1) that is due to the softness of the silicone-
rubber material[21].

Figure 1: (left) Top view and (right) side view of the experimental setup. The
metafluid is confined in a prism-shaped box, with an angle
$\theta_{\rm{fluid}}=+2$°, that is deposited directly on a large broad-band
ultrasonic immersion transducer (150 mm x 40 mm). The acoustic waves propagate
from the bottom to the top of the water tank, along the vertical $z$-axis. The
acoustic field refracted from the prism is scanned in the $x$-$z$ plane (33 mm
x 33 mm) with a small ultrasonic probe.

## Results

When excited with a short electrical pulse, the large ultrasonic transducer
generates a short broad-band acoustic pulse propagating in the confined fluid
along the $z$-axis only. At the interface between the confined fluid and
water, the acoustic pulse is then refracted in the surrounding water with an
angle of refraction $\theta_{\rm{water}}$ shown in Fig. 1. By measuring
$\theta_{\rm{water}}$, the acoustic index $n_{\rm{fluid}}$ of the fluid
confined in the prism can be easily deduced from the Snell’s law as following:

$n_{\rm{fluid}}=n_{\rm{water}}\frac{\rm{sin}(\theta_{\rm{water}})}{\rm{sin}(\theta_{\rm{fluid}})}$
(1)

with $n_{\rm{water}}=1$ since the acoustic index $n$ ($=c_{0}/c$) of a
material with the phase velocity $c$ is usually defined relatively to water
($c_{0}=c_{\textrm{water}}$) for underwater acoustics. In these experiments,
the acoustic index $n_{\rm{fluid}}$ cannot be directly retrieved from the
refracted transmitted temporal signals measured in the $x$-$z$ plane as
depicted in Fig. 1, since a short broad-band acoustic pulse has been used in
these pulsed ultrasonic experiments. In a such broad frequency range (from 15
kHz to 600 kHz), the acoustic index of a concentrated metafluid is expected to
exhibit strong variations ranging from high positive values to negative
ones[15]. Due to potential strong dispersion effects, we focus here on a
harmonic spectral analysis of the refracted beams by performing Fourier
transforms over time and space. First, we performed time-domain Fourier
transforms of all the signals acquired at each position over a 2D spatial grid
(Fig. 2.a.) in order to obtain the spatial-field at each frequency component
(Fig. 2.b.). Then a 2D spatial Fourier transform is applied at each frequency
for getting the wavenumber spectrum of the harmonic beam (Fig. 2.c.). Since
the beams refracted in water are quasi-plane waves, it is straightforward to
extract accurately their direction of propagation from the (real-valued)
spatial Fourier components $k^{\textrm{r}}_{\textrm{x}}$ and
$k^{\textrm{r}}_{\textrm{z}}$ at the peak amplitude of the wavenumber
spectrum.

Figure 2: (a) Snapshot of the measured pressure field refracted from the prism
filled with concentrated metafluid ($\Phi=17\%$). (b) Corresponding angular
phase, shown here at 175 kHz, obtained from the time-domain Fourier transforms
performed for each position over the square grid shown in (a). (c) 2D spatial
Fourier transform of the scanned field pattern at 175 kHz for the measurement
of the (real-valued) spatial Fourier components $k^{\textrm{r}}_{\textrm{x}}$
and $k^{\textrm{r}}_{\textrm{z}}$ of the refracted beam in water. Figure 3: 2D
Fourier transforms of the scanned field pattern performed at different
frequencies: 40, 130, 160, 200, 250 kHz (from the top to bottom). The prism is
filled with either water (left), diluted metafluid (center) or concentrated
metafluid (right). On all maps, the white vertical axis ($k_{\textrm{x}}=0$)
corresponds to the case in which
$\theta_{\textrm{water}}=\theta_{\textrm{fluid}}$ meaning that
$n_{\textrm{fluid}}=+1$, whereas the tilted white axis
($\arctan(k_{\textrm{x}}/k_{\textrm{z}})=-\theta_{\rm{fluid}}$) corresponds to
the case in which $\theta_{\textrm{water}}=0$ meaning that
$n_{\textrm{fluid}}=0$. The Fourier transforms imply a plane wave in the
positive $z$ direction.

From the acoustic field map shown at a given time in Fig. 2.a, we can get the
corresponding map of the angular phase for each frequency component of the
refracted beam. As an example, Fig. 2.b shows the angular phase of the
refracted beam at $f=175$ kHz revealing an angle of refraction of a few
degrees from the initially vertical propagation direction. Such a deviation
can be estimated by measuring the tilted angle of these refracted wavefronts
in the $x$-$z$ plane but this direct spatial measurement may suffer from
fluctuations observed on the wavefronts. As an alternative, we performed 2D
Fourier transforms of the scanned field patterns at different frequencies to
infer the $k^{\textrm{r}}_{\textrm{x}}$ and $k^{\textrm{r}}_{\textrm{z}}$
coordinates of the wave vector $\textbf{k}_{\textrm{water}}$ in the
$k_{\textrm{x}}$-$k_{\textrm{z}}$ plane. The values of these two components
are given by the location of the maximum spot intensity shown in Fig. 2.c.
Then, the angle of refraction $\theta_{\rm{water}}$ can be easily deduced for
each frequency as following:

$\theta_{\rm{water}}=\arctan(\frac{k^{\textrm{r}}_{\textrm{x}}}{k^{\textrm{r}}_{\textrm{z}}})+\theta_{\rm{fluid}}$
(2)

In these experiments, the coordinate $k^{\textrm{r}}_{\textrm{z}}$ is
necessarily positive since the acoustic waves propagate from the bottom to the
top of water tank along the $z$-axis in any case. However, the coordinate
$k^{\textrm{r}}_{\textrm{x}}$ may be either positive or negative depending on
the direction of the refracted beam. The acoustic index $n_{\rm{fluid}}$ of
the fluid confined in the prism can be deduced from the values of the angle
$\theta_{\rm{water}}$ by using Eq. 1. When the prism is filled with water
(Fig. 3, left), the spot goes along the vertical $k_{z}$-axis for which
$k^{\textrm{r}}_{\textrm{x}}=0$ (white vertical line) as the frequency
increases leading to $\theta_{\rm{water}}=\theta_{\rm{fluid}}$ according to
Eq. 2. Therefore, no refraction occurs in that case which may be expected
since the material is the same (water) on both sides of the interface. Note
that if the spot had gone along the white titled line shown in Figs. 3, for
which
$\arctan(k^{\textrm{r}}_{\textrm{x}}/k^{\textrm{r}}_{\textrm{z}})=-\theta_{\textrm{fluid}}$,
this would lead to $\theta_{\textrm{water}}=0$, corresponding thus to a zero-
index fluid confined in the prism.

When the prism is filled with the diluted metafluid ($\Phi=3\%$, $d=750$
$\mu$m with a size dispersion of 30%), the spot oscillates around the white
vertical axis as the frequency increases (Fig. 3, center). The values of
$n_{\textrm{fluid}}$ extracted from Eqs.1 and 2 show that the acoustic index
of this diluted metafluid varies from high values at low frequencies
($n_{\textrm{fluid}}=+4$ at 40 kHz) to low values at intermediate frequencies
($n_{\textrm{fluid}}=+0.5$ around 130 kHz) before getting closer to +1 at high
frequencies as shown in Fig. 4. This strong dispersion is due to the micro-
bead resonances that occur around 150 kHz for that size of particles. Far away
from these low-frequency acoustic resonances, the acoustic index of the
metafluid is similar to that of the aqueous matrix[15].

We also produced a concentrated metafluid ($\Phi=17\%$, $d=700$ $\mu$m with a
size dispersion of 10%) in which negative refraction occurs for a certain
range of frequencies. Actually, Fig. 3 (right) shows that the spot can go
beyond the white tilted axis (corresponding to $n_{\textrm{fluid}}=0$) as
observed at 160 kHz. The experimental results are not shown at 40 kHz because
of the very high attenuation that is due the strong (monopolar) low-frequency
resonances of the micro-beads [15]. But the acoustic index is shown to be
negative over a broad frequency range as shown in Fig. 4. Note that this
’negative band’ is slightly shifted to higher frequencies compared to the
diluted metafluid because of the smaller size of the micro-beads of this
concentrated metafluid [22].

Figure 4: Acoustic index $n_{\textrm{fluid}}$ for different fluids confined in
the prism-shaped box and extracted as a function of frequency from the
$k$-maps. Dashed lines refer to calculations performed in the framework of
multiple scattering theory (see the text).

Finally, we compared our acoustical measurements to theoretical predictions
produced through multiple-scattering modeling, revealing good qualitative
agreement as shown in Fig. 4. The theoretical acoustic index of the metafluids
($\Phi=3\%$ and $17\%$) are obtained from their effective wavenumbers
calculated using the Waterman-Truell formula[23]. The values of the material
parameters for the soft porous silicone rubber that we used for these
calculations, were $c_{L}$ = 120 m.s-1, $\alpha_{L}$= 20 Np.m-1.MHz-2 (the
phase velocity and attenuation coefficient for longitudinal waves), $c_{T}$ =
40 m.s-1, $\alpha_{T}$= 200 Np.m-1.MHz-2 (the phase velocity and attenuation
coefficient for shear waves) and $\rho_{1}$ = 760 kg.m-3. The water-based gel
matrix has the same properties as water ($\rho_{0}$ = 1000 kg.m-3, $c_{0}$ =
1490 m.s-1) since this host Bingham fluid is essentially made of water and
very small amounts of polymer (Carbopol) to prevent the creaming of the porous
micro-beads.

## Conclusion

In summary, we have reported an experimental demonstration of negative
refraction in a 3D acoustic metamaterial. The experiments have been conducted
at ultrasonic frequencies with a locally resonant metafluid composed of soft
porous silicone-rubber micro-beads whose concentration must be high enough so
that refraction be negative. However, the achievement of acoustic devices
based on negative refraction such as perfect lenses envisioned by Pendry[24]
seems unattainable with our metafluids because of their large attenuation that
is mainly due to the strong resonant scattering[15]. Alternatively to bulky
and lossy 3D metamaterials, acoustic metasurfaces[25] might be much more
appropriate to manipulate acoustic wavefronts as recently demonstrated with
soft gradient-index metasurfaces[26].

## Methods

In this study, we used a large ultrasonic transducer (150 mm x 40 mm), with a
central frequency of 150 kHz, to ensure the generation of quasi-pure plane
waves in the water tank (propagating along the vertical axis from the bottom
to the top of the water tank). This broad-band transducer was excited with a
short electrical pulse generated by a pulser/receiver (Olympus, 5077PR) that
was also used to amplify the electric signal recorded by the 1-inch-diameter
receiving transducer (Olympus V301) before its acquisition on a computer via a
waveform digitizer card (AlazarTech, ATS460). The angle of the prism has been
chosen here as low as possible ($\theta_{\rm{fluid}}=+2$°) in order to
guarantee that the refracted-beam amplitude does not vary too much along the
x-axis in spite of the slightly increasing depth of the prism. For the
experiments conducted with the prism filled with water (Fig. 3, left) or with
the diluted metafluid (Fig. 3, centre), the width of the scan area along the
x-axis was 60 mm. For the concentrated metafluid (Fig. 3, right), this width
was reduced to 30 mm inducing a slight spread of the Fourier Transform
distribution along the $k_{x}$-axis. The acoustical fields refracted in water
were scanned on grids with a step of 1 mm that is 3 times smaller than the
acoustic wavelength in water at 500 kHz (= 3 mm). Note that we also used zero-
padding techniques before doing Fast Fourier Transforms in space to improve
the resolution in the $k$-maps shown in Fig. 3.

## References

  * [1] Liu, Z. _et al._ Locally resonant sonic materials.  _Science_ 289, 1734–1736 (2000).
  * [2] Ma, G. & Sheng, P.  Acoustic metamaterials: from local resonances to broad horizons.  _Sci. Adv._ 2, e1501595 (2016).
  * [3] Cummer, S. A., Christensen, J. & Alù, A.  Controlling sound with acoustic metamaterials.  _Nat. Rev. Mater._ 1, 16001 EP – (2016).
  * [4] Li, J. & Chan, C. T.  Double-negative acoustic metamaterial.  _Phys. Rev. E_ 70, 055602 (2004).
  * [5] Lee, S. H., Park, C. M., Seo, Y. M., Wang, Z. G. & Kim, C. K.  Composite acoustic medium with simultaneously negative density and modulus.  _Phys. Rev. Lett._ 104, 054301 (2010).
  * [6] Liang, Z., Willatzen, M., Li, J. & Christensen, J.  Tunable acoustic double negativity metamaterial.  _Sci. Rep._ 2, 859 (2012).
  * [7] Yang, M., Ma, G., Yang, Z. & Sheng, P.  Coupled membranes with doubly negative mass density and bulk modulus.  _Phys. Rev. Lett._ 110, 134301 (2013).
  * [8] Maurya, S. K., Pandey, A., Shukla, S. & Saxena, S.  Double negativity in 3d space coiling metamaterials.  _Sci. Rep._ 6, 33683 (2016).
  * [9] Lanoy, M. _et al._ Acoustic double negativity induced by position correlations within a disordered set of monopolar resonators.  _Phys. Rev. B_ 96, 220201 (2017).
  * [10] Zhou, Y., Fang, X., Li, D., Hao, T. & Li, Y.  Acoustic multiband double negativity from coupled single-negative resonators.  _Phys. Rev. Applied_ 10, 044006 (2018).
  * [11] Wang, W., Bonello, B., Djafari-Rouhani, B., Pennec, Y. & Zhao, J.  Double-negative pillared elastic metamaterial.  _Phys. Rev. Applied_ 10, 064011 (2018).
  * [12] Dong, H.-W., Zhao, S.-D., Wang, Y.-S., Cheng, L. & Zhang, C.  Robust 2d/3d multi-polar acoustic metamaterials with broadband double negativity.  _J. Mech. Phys. Solids_ 137, 103889 (2020).
  * [13] Li, J., Fung, K. H., Liu, Z. Y., Sheng, P. & Chan, C. T.  _Physics of Negative Refraction and Negative Index Materials_.  edited by C. M. Krowne and Y. Zong (Springer, Berlin, 2007).
  * [14] Brunet, T., Poncelet, O. & Aristégui, C.  Negative-index metamaterials: is double negativity a real issue for dissipative media?  _EPJ Appl. Metamat._ 2, 3 (2015).
  * [15] Brunet, T. _et al._ Soft 3D acoustic metamaterial with negative index.  _Nat. Mater._ 14, 384–388 (2015).
  * [16] Kaina, N., Lemoult, F., Fink, M. & Lerosey, G.  Negative refractive index and acoustic superlens from multiple scattering in single negative metamaterials.  _Nature_ 525, 77–81 (2015).
  * [17] Brunet, T., Leng, J. & Mondain-Monval, O.  Soft acoustic metamaterials.  _Science_ 342, 323–324 (2013).
  * [18] Raffy, S., Mascaro, B., Brunet, T., Mondain-Monval, O. & Leng, J.  A soft 3d acoustic metafluid with dual-band negative refractive index.  _Adv. Mater._ 28, 1760–1764 (2016).
  * [19] Popa, B.-I. & Cummer, S. A.  Water-based metamaterials: Negative refraction of sound.  _Nat. Mater._ 14, 363–364 (2015).
  * [20] Kovalenko, A., Zimny, K., Mascaro, B., Brunet, T. & Mondain-Monval, O.  Tailoring of the porous structure of soft emulsion-templated polymer materials.  _Soft Matter_ 12, 5154–5163 (2016).
  * [21] Ba, A., Kovalenko, A., Aristégui, C., Mondain-Monval, O. & Brunet, T.  Soft porous silicone rubbers with ultra-low sound speeds in acoustic metamaterials.  _Sci. Rep._ 7, 40106 EP – (2017).
  * [22] Brunet, T. _et al._ Sharp acoustic multipolar-resonances in highly monodisperse emulsions.  _Appl. Phys. Lett._ 101, 011913 (2012).
  * [23] Waterman, P. C. & Truell, R.  Multiple scattering of waves.  _J. Math. Phys._ 2, 512–537 (1961).
  * [24] Pendry, J. B.  Negative refraction makes a perfect lens.  _Phys. Rev. Lett._ 85, 3966–3969 (2000).
  * [25] Assouar, B. _et al._ Acoustic metasurfaces.  _Nat. Rev. Mater._ 3, 460–472 (2018).
  * [26] Jin, Y., Kumar, R., Poncelet, O., Mondain-Monval, O. & Brunet, T.  Flat acoustics with soft gradient-index metasurfaces.  _Nat. Commun._ 10, 143 (2019).

## Acknowledgements

This work was partially funded and performed within the framework of the Labex
AMADEUS ANR-10-LABEX-0042-AMADEUS with the help of the French state Initiative
d’Excellence IdEx ANR-10-IDEX-003-02 and project BRENNUS ANR-15-CE08-0024 (ANR
and FRAE funds).

## Author contributions statement

T.B. supervised the project, A.K. produced the soft porous silicone-rubber
micro-beads to achieve the metafluids under the guidance of O.M.-M., B.T. and
T.B. conducted the underwater experiments, B.T., T.B., C.A., and O.P analyzed
the results. All authors reviewed the manuscript.

## Additional information

Competing financial interests The authors declare no competing financial
interests.