<figure><img src="image_1.png" /></figure>

HAL
open science

# Remote sensing and hydrogeophysics give a new impetus to integrated hydrological models: A review

M.W. Lubczynski, M. Leblanc, O. Batelaan

► To cite this version:

M.W. Lubczynski, M. Leblanc, O. Batelaan. Remote sensing and hydrogeophysics give a new impetus to integrated hydrological models: A review. Journal of Hydrology, 2024, 633, pp.130901. 10.1016/j.jhydrol.2024.130901. hal-04573400

HAL Id: hal-04573400

https://hal.inrae.fr/hal-04573400

Submitted on 13 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Remote sensing and hydrogeophysics give a new impetus to integrated hydrological models: A review

M.W. Lubczynski, M. Leblanc, O. Batelaan

► To cite this version:

M.W. Lubczynski, M. Leblanc, O. Batelaan. Remote sensing and hydrogeophysics give a new impetus to integrated hydrological models: A review. Journal of Hydrology, 2024, 633, pp.130901. 10.1016/j.jhydrol.2024.130901. hal-04573400

HAL Id: hal-04573400

https://hal.inrae.fr/hal-04573400v1

Submitted on 13 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Review papers

# Remote sensing and hydrogeophysics give a new impetus to integrated hydrological models: A review

M.W. Lubczynski<sup>a, *</sup>, M. Leblanc<sup>b, c</sup>, O. Batelaan<sup>d</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

This manuscript was handled by Corrado Corradini, Editor-in-Chief, with the assistance of Carla Saltalippi, Associate Editor

**Keywords:**

Remote sensing and Hydrogeophysics

Climate forcings

Parameters

Boundary conditions

Observations to constrain simulated state

variables

Integrated hydrological models

## ABSTRACT

Integrated Hydrological Models (IHMs) dynamically couple surface and groundwater processes across the unsaturated zone domain. IHMs are data intensive and computationally demanding but can provide physically realistic output, particularly if sufficient input data of high quality is available. In-situ observations often have a small footprint and are time and cost-demanding. Satellite remote sensing observations, with their long time series archives and spatially semi-continuous gridded format, as well as hydrogeophysical observations with their flexible, 'on-demand' high-resolution data coverage, perfectly complement in-situ observations. We review the contribution of various satellite remote sensing products for IHM: (1) climate forcings, (2) parameters, (3) boundary conditions and (4) observations for constraining model calibration and data assimilation. Our review of hydrogeophysics focuses on the four mentioned IHM contributions, but we analyze them per data acquisition platform, i.e., surface, drone-borne and airborne hydrogeophysics. Finally, the review includes a discussion on the optimal use of satellite remote sensing and hydrogeophysical data in IHMs, as well as a vision for further improvements of data-driven, integrated hydrological modelling.

## 1. Introduction

Modern water management requires spatially distributed numerical hydrological models as tools to integrate data for analyses and forecasting scenarios. Traditionally, hydrologists have focused on modelling specific hydrological domains, i.e., the saturated zone, unsaturated zone or land surface-atmosphere (including hydrological processes in open water bodies). Consequently, different modelling expertises were developed to address these hydrological domains, i.e., land surface modelling (including 2D watershed modelling, also described as rainfall-runoff models) and subsurface modelling, the latter separately dealing with either unsaturated zone or saturated zone domains. The models addressing only one hydrological domain are referred to hereafter as standalone models. Various simplifying assumptions are made in each model type, generating conceptual or structural uncertainty (Enemark et al., 2019).

A land surface model uses quantitative methods to simulate the exchange of water and energy fluxes at the Earth surface-atmosphere interface. Some of these models, such as the rainfall-runoff models,

can simulate streamflow generation and routing based on meteorological data and land surface characteristics. The main output is generally the streamflow hydrograph at the outlet of a simulated watershed, which is compared with the measured hydrograph during model calibration. Also, evapotranspiration (ET) is an output of land surface models, so it can be compared or calibrated with ground ET measurements, e.g., with data from eddy covariance towers or ET observations from space. Like other models, land surface models have inherent uncertainties due to system heterogeneity. However, their largest uncertainty is at their bottom boundary. Most do not include or simplify subsurface flow, neglecting the subsurface lateral flow (Zhao et al., 2021) and flow exchanges between land surface and aquifers. In contrast, in groundwater models, the most uncertain boundary is the top boundary, i.e., the spatiotemporally variable net recharge at the interface between the unsaturated and the saturated zone. Unsaturated zone models are the most uncertain because their top and bottom boundaries are uncertain, and also, their variably-saturated domain requires solutions of highly non-linear flow equations.

<footnote><sup>a</sup> Department of Water Resources, Faculty of Geo-Information Science and Earth Observation (ITC), University of Twente, PO Box 217, 7500 AE Enschede, the Netherlands</footnote>

<footnote><sup>b</sup> INRAE, Avignon University, UMR EMMAH, F-84000 Avignon, France</footnote>

<footnote><sup>c</sup> IRD, UMR G-EAU, University Mohammed V, Geosciences Laboratory, Rabat, Morocco</footnote>

<footnote><sup>d</sup> National Centre for Groundwater Research and Training (NCGRT), College of Science and Engineering, Flinders University, Adelaide, SA 5001, Australia</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: m.w.lubczynski@utwente.nl, lubczynskimaciek@gmail.com (M.W. Lubczynski).</footnote>
Review papers

# Remote sensing and hydrogeophysics give a new impetus to integrated hydrological models: A review

M.W. Lubczynski<sup>a, *</sup>, M. Leblanc<sup>b, c</sup>, O. Batelaan<sup>d</sup>

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

This manuscript was handled by Corrado Corradini, Editor-in-Chief, with the assistance of Carla Saltalippi, Associate Editor

**Keywords:**

Remote sensing and Hydrogeophysics

Climate forcings

Parameters

Boundary conditions

Observations to constrain simulated state

variables

Integrated hydrological models

## ABSTRACT

Integrated Hydrological Models (IHMs) dynamically couple surface and groundwater processes across the unsaturated zone domain. IHMs are data intensive and computationally demanding but can provide physically realistic output, particularly if sufficient input data of high quality is available. In-situ observations often have a small footprint and are time and cost-demanding. Satellite remote sensing observations, with their long time series archives and spatially semi-continuous gridded format, as well as hydrogeophysical observations with their flexible, 'on-demand' high-resolution data coverage, perfectly complement in-situ observations. We review the contribution of various satellite remote sensing products for IHM: (1) climate forcings, (2) parameters, (3) boundary conditions and (4) observations for constraining model calibration and data assimilation. Our review of hydrogeophysics focuses on the four mentioned IHM contributions, but we analyze them per data acquisition platform, i.e., surface, drone-borne and airborne hydrogeophysics. Finally, the review includes a discussion on the optimal use of satellite remote sensing and hydrogeophysical data in IHMs, as well as a vision for further improvements of data-driven, integrated hydrological modelling.

## 1. Introduction

Modern water management requires spatially distributed numerical hydrological models as tools to integrate data for analyses and forecasting scenarios. Traditionally, hydrologists have focused on modelling specific hydrological domains, i.e., the saturated zone, unsaturated zone or land surface-atmosphere (including hydrological processes in open water bodies). Consequently, different modelling expertises were developed to address these hydrological domains, i.e., land surface modelling (including 2D watershed modelling, also described as rainfall-runoff models) and subsurface modelling, the latter separately dealing with either unsaturated zone or saturated zone domains. The models addressing only one hydrological domain are referred to hereafter as standalone models. Various simplifying assumptions are made in each model type, generating conceptual or structural uncertainty (Enemark et al., 2019).

A land surface model uses quantitative methods to simulate the exchange of water and energy fluxes at the Earth surface-atmosphere interface. Some of these models, such as the rainfall-runoff models,

can simulate streamflow generation and routing based on meteorological data and land surface characteristics. The main output is generally the streamflow hydrograph at the outlet of a simulated watershed, which is compared with the measured hydrograph during model calibration. Also, evapotranspiration (ET) is an output of land surface models, so it can be compared or calibrated with ground ET measurements, e.g., with data from eddy covariance towers or ET observations from space. Like other models, land surface models have inherent uncertainties due to system heterogeneity. However, their largest uncertainty is at their bottom boundary. Most do not include or simplify subsurface flow, neglecting the subsurface lateral flow (Zhao et al., 2021) and flow exchanges between land surface and aquifers. In contrast, in groundwater models, the most uncertain boundary is the top boundary, i.e., the spatiotemporally variable net recharge at the interface between the unsaturated and the saturated zone. Unsaturated zone models are the most uncertain because their top and bottom boundaries are uncertain, and also, their variably-saturated domain requires solutions of highly non-linear flow equations.

<footnote><sup>a</sup> Department of Water Resources, Faculty of Geo-Information Science and Earth Observation (ITC), University of Twente, PO Box 217, 7500 AE Enschede, the Netherlands</footnote>

<footnote><sup>b</sup> INRAE, Avignon University, UMR EMMAH, F-84000 Avignon, France</footnote>

<footnote><sup>c</sup> IRD, UMR G-EAU, University Mohammed V, Geosciences Laboratory, Rabat, Morocco</footnote>

<footnote><sup>d</sup> National Centre for Groundwater Research and Training (NCGRT), College of Science and Engineering, Flinders University, Adelaide, SA 5001, Australia</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: m.w.lubczynski@utwente.nl, lubczynskimaciek@gmail.com (M.W. Lubczynski).</footnote>
In reality, the saturated, unsaturated and surface flow domains represent one dynamically interconnected hydrological system. Integrated hydrological models (IHMs) combine these three hydrological domains in different ways depending on the IHM type into one dynamic system. The IHMs are forced by external climate forcings (e.g. rainfall), restricted by boundary conditions and regulated by parameters of all three domains, which are adjusted in a model calibration by minimizing the difference between calibrated state variables and observations. The term 'calibrated state variable' or just 'state variable' was assigned after Ahlfeld et al. (2011) and represents the modelled state of the hydrological system at a particular time and location, for example, hydraulic head, temperature, soil moisture, storage or concentration, but also evapotranspiration flux or river flow. As such, next to system conceptualisation, the reliability of IHMs depends on the quality of a wide range of data, spatiotemporal data coverage, the number of different types of observations applied for calibration of state variables, and the goodness of model calibration.

The acquisition of in-situ point data and their use for setting up and calibrating models is important but time-consuming, expensive and usually spatially-wise insufficient. Point data typically have small footprints, i.e. the surface area or subsurface volume from which the measured signal arises. So, they have limited capacity to represent spatiotemporal variability. Even a large number of in-situ point data distributed over a study area can be insufficient because interpolation usually cannot reproduce its true spatiotemporal variability. In that respect, climate forcings of IHMs are the most vulnerable to interpolation errors, because they are the most spatiotemporally variable (especially precipitation), but have to be assigned for each stress period of a model. Remote sensing (RS) and hydrogeophysics, when bias-corrected using in-situ point data and applying appropriate geostatistical methods, can better handle the spatiotemporal variability of model input data than interpolation alone. Therefore, RS and hydrogeophysics are increasingly being used as sources of input data for distributed models, including IHMs. This manuscript aims to review current and potential RS and hydrogeophysics contributions to IHMs to increase their reliability and stimulate wider use of RS and hydrogeophysics data in IHMs.

Remote sensing and hydrogeophysics have the same main principle, i.e., they acquire information about the target from a distance without direct contact with it. This is why the scientific literature has terminology-related unclarity on where to put the boundary between RS and hydrogeophysics. Moreover, there are even opinions that geophysics is part of remote sensing or that remote sensing is part of geophysics. In this study, we decided to compromise on the terminology. Traditionally, sensing the ground surface and subsurface from satellites has always been referred to as remote sensing, while sensing the subsurface from the ground has always been referred to as surface geophysics or surface hydrogeophysics if targeting water. This tradition is respected in this review, and the terms satellite remote sensing or just remote sensing refer exclusively to data acquisition from satellites, while the data acquisition from the ground surface is referred to as surface hydrogeophysics. The remaining data acquisition type from above, but from near the ground surface, is where the largest terminological controversy is present. Unlike most satellite data acquisition programs, these data acquisitions are typically made 'on demand' (as in surface geophysics) and also from a much lower height. Hence, they cover much smaller footprints than satellites but at higher spatial resolution. As such we classify them as hydrogeophysics, being divided into airborne hydrogeophysics, where data acquisition occurs from pilot-operated flying machines (from helicopters or small airplanes) and drone-borne hydrogeophysics, where data acquisition occurs from drones, remotely operated from the ground surface. Note, only non-invasive hydrogeophysical methods are discussed hereafter.

## 2. Integrated hydrological models

In the last decades, many physically based modelling approaches have been developed addressing multiple processes interacting over a range of scales and across different hydrological domains (Paniconi and Putti, 2015), giving rise to integrated hydrological models. Each IHM has its own way of adapting the governing equations for the surface and the subsurface flow domains, the coupling strategy and the solution technique (Haque et al., 2021; Maxwell et al., 2014; Ntona et al., 2022; Yang et al., 2021). The most important characteristic of IHMs is that they dynamically couple the surface water and groundwater flow domains across the unsaturated zone, meaning that the solution of several flow domains is computed within the same single model simulation.

Following Daoud et al. (2022), IHMs can be divided according to the complexity of integration of the surface water and groundwater flow domains into: (1) fully-coupled IHMs; and (2) simplified IHMs. In the fully-coupled IHMs, the physically-based governing equations of all modelling domains are solved simultaneously. Examples of fully-coupled IHM codes are CATHY (Camporese et al., 2010), HydroGeoSphere (Therrien and Sudicky, 1996), MODHSM (Panday and Huyakorn, 2004), and PARFLOW (Kollet and Maxwell, 2006). Those fully-coupled IHMs provide a high level of process detail, but are also costly in terms of modelling time and computational requirements. The simplified IHMs, also dynamically couple surface water and groundwater, but simplify governing equations of one or more flow domains, achieving a computationally more efficient model than a fully-coupled IHM. Examples of simplified IHM codes are various MODFLOW codes extended with the UZF Package, where vertical flow across the unsaturated zone is simulated by the kinematic wave approximation (Niswonger et al., 2006) of the 1D Richards' equation. The IHMs with such simplification require fewer input parameters than those including a solution of the Richards' equation and avoid the highly non-linear and computationally intensive solution of the Richards' equation. Thanks to their efficiency, the simplified IHMs are particularly suitable for catchment or coarser-scale model applications, because the potential errors associated with the kinematic wave approximation are small relative to the errors resulting from scaling effects and from a reduced set of parameters representative of a highly complex system (Morway et al., 2012).

The simplified IHMs can be either: (a) multi-environment; or (b) single-environment. A multi-environment IHM involves the dynamic coupling of two or more codes, such as the coupling of MODFLOW with PRMS in the GSFLOW (Markstrom et al., 2008). A single-environment IHM code, e.g. MODFLOW 6 (Hughes et al., 2017), internally simulates flow across the unsaturated zone by using so-called advanced packages (Langevin et al., 2017), which dynamically integrate the surface flow domain with the groundwater flow domain across the unsaturated zone domain. It can also integrate different model components through the Water Mover (MVR) Package, which allows the transfer of water fluxes between different hydrological domains (Daoud et al., 2022).

Even if well-calibrated using ample, high-quality input data, all hydrological models still have inherent uncertainty. The data assimilation technique, commonly used in weather forecasting, where model states are continuously updated with real-time data, is increasingly used in IHMs. It allows not only to adapt a model to real-time states but also to improve model parameterisation with new data, reducing model uncertainty (Camporese and Girotto, 2022; De Lannoy et al., 2022; Doherty and Moore, 2020; Liu et al., 2012). Considering groundwater IHM studies, so far, the data assimilation has gained mainly interest through research applications, applying complex multi-environment, Richards' equation-based models (He et al., 2019; Zhang et al., 2018). However, recently, more efficient, practical engineering solutions have also become available (Doherty and Moore, 2020), such as the open-source data assimilation tool PESTPP-DA, which is part of the PEST++ framework (White et al., 2020), run under the pyEMU environment (White et al., 2016).

## 3. Application of remote sensing in integrated hydrological models

RS can contribute to IHMs by assigning: i) climate forcings, also
In reality, the saturated, unsaturated and surface flow domains represent one dynamically interconnected hydrological system. Integrated hydrological models (IHMs) combine these three hydrological domains in different ways depending on the IHM type into one dynamic system. The IHMs are forced by external climate forcings (e.g. rainfall), restricted by boundary conditions and regulated by parameters of all three domains, which are adjusted in a model calibration by minimizing the difference between calibrated state variables and observations. The term 'calibrated state variable' or just 'state variable' was assigned after Ahlfeld et al. (2011) and represents the modelled state of the hydrological system at a particular time and location, for example, hydraulic head, temperature, soil moisture, storage or concentration, but also evapotranspiration flux or river flow. As such, next to system conceptualisation, the reliability of IHMs depends on the quality of a wide range of data, spatiotemporal data coverage, the number of different types of observations applied for calibration of state variables, and the goodness of model calibration.

The acquisition of in-situ point data and their use for setting up and calibrating models is important but time-consuming, expensive and usually spatially-wise insufficient. Point data typically have small footprints, i.e. the surface area or subsurface volume from which the measured signal arises. So, they have limited capacity to represent spatiotemporal variability. Even a large number of in-situ point data distributed over a study area can be insufficient because interpolation usually cannot reproduce its true spatiotemporal variability. In that respect, climate forcings of IHMs are the most vulnerable to interpolation errors, because they are the most spatiotemporally variable (especially precipitation), but have to be assigned for each stress period of a model. Remote sensing (RS) and hydrogeophysics, when bias-corrected using in-situ point data and applying appropriate geostatistical methods, can better handle the spatiotemporal variability of model input data than interpolation alone. Therefore, RS and hydrogeophysics are increasingly being used as sources of input data for distributed models, including IHMs. This manuscript aims to review current and potential RS and hydrogeophysics contributions to IHMs to increase their reliability and stimulate wider use of RS and hydrogeophysics data in IHMs.

Remote sensing and hydrogeophysics have the same main principle, i.e., they acquire information about the target from a distance without direct contact with it. This is why the scientific literature has terminology-related unclarity on where to put the boundary between RS and hydrogeophysics. Moreover, there are even opinions that geophysics is part of remote sensing or that remote sensing is part of geophysics. In this study, we decided to compromise on the terminology. Traditionally, sensing the ground surface and subsurface from satellites has always been referred to as remote sensing, while sensing the subsurface from the ground has always been referred to as surface geophysics or surface hydrogeophysics if targeting water. This tradition is respected in this review, and the terms satellite remote sensing or just remote sensing refer exclusively to data acquisition from satellites, while the data acquisition from the ground surface is referred to as surface hydrogeophysics. The remaining data acquisition type from above, but from near the ground surface, is where the largest terminological controversy is present. Unlike most satellite data acquisition programs, these data acquisitions are typically made 'on demand' (as in surface geophysics) and also from a much lower height. Hence, they cover much smaller footprints than satellites but at higher spatial resolution. As such we classify them as hydrogeophysics, being divided into airborne hydrogeophysics, where data acquisition occurs from pilot-operated flying machines (from helicopters or small airplanes) and drone-borne hydrogeophysics, where data acquisition occurs from drones, remotely operated from the ground surface. Note, only non-invasive hydrogeophysical methods are discussed hereafter.

## 2. Integrated hydrological models

In the last decades, many physically based modelling approaches have been developed addressing multiple processes interacting over a range of scales and across different hydrological domains (Paniconi and Putti, 2015), giving rise to integrated hydrological models. Each IHM has its own way of adapting the governing equations for the surface and the subsurface flow domains, the coupling strategy and the solution technique (Haque et al., 2021; Maxwell et al., 2014; Ntona et al., 2022; Yang et al., 2021). The most important characteristic of IHMs is that they dynamically couple the surface water and groundwater flow domains across the unsaturated zone, meaning that the solution of several flow domains is computed within the same single model simulation.

Following Daoud et al. (2022), IHMs can be divided according to the complexity of integration of the surface water and groundwater flow domains into: (1) fully-coupled IHMs; and (2) simplified IHMs. In the fully-coupled IHMs, the physically-based governing equations of all modelling domains are solved simultaneously. Examples of fully-coupled IHM codes are CATHY (Camporese et al., 2010), HydroGeoSphere (Therrien and Sudicky, 1996), MODHSM (Panday and Huyakorn, 2004), and PARFLOW (Kollet and Maxwell, 2006). Those fully-coupled IHMs provide a high level of process detail, but are also costly in terms of modelling time and computational requirements. The simplified IHMs, also dynamically couple surface water and groundwater, but simplify governing equations of one or more flow domains, achieving a computationally more efficient model than a fully-coupled IHM. Examples of simplified IHM codes are various MODFLOW codes extended with the UZF Package, where vertical flow across the unsaturated zone is simulated by the kinematic wave approximation (Niswonger et al., 2006) of the 1D Richards' equation. The IHMs with such simplification require fewer input parameters than those including a solution of the Richards' equation and avoid the highly non-linear and computationally intensive solution of the Richards' equation. Thanks to their efficiency, the simplified IHMs are particularly suitable for catchment or coarser-scale model applications, because the potential errors associated with the kinematic wave approximation are small relative to the errors resulting from scaling effects and from a reduced set of parameters representative of a highly complex system (Morway et al., 2012).

The simplified IHMs can be either: (a) multi-environment; or (b) single-environment. A multi-environment IHM involves the dynamic coupling of two or more codes, such as the coupling of MODFLOW with PRMS in the GSFLOW (Markstrom et al., 2008). A single-environment IHM code, e.g. MODFLOW 6 (Hughes et al., 2017), internally simulates flow across the unsaturated zone by using so-called advanced packages (Langevin et al., 2017), which dynamically integrate the surface flow domain with the groundwater flow domain across the unsaturated zone domain. It can also integrate different model components through the Water Mover (MVR) Package, which allows the transfer of water fluxes between different hydrological domains (Daoud et al., 2022).

Even if well-calibrated using ample, high-quality input data, all hydrological models still have inherent uncertainty. The data assimilation technique, commonly used in weather forecasting, where model states are continuously updated with real-time data, is increasingly used in IHMs. It allows not only to adapt a model to real-time states but also to improve model parameterisation with new data, reducing model uncertainty (Camporese and Girotto, 2022; De Lannoy et al., 2022; Doherty and Moore, 2020; Liu et al., 2012). Considering groundwater IHM studies, so far, the data assimilation has gained mainly interest through research applications, applying complex multi-environment, Richards' equation-based models (He et al., 2019; Zhang et al., 2018). However, recently, more efficient, practical engineering solutions have also become available (Doherty and Moore, 2020), such as the open-source data assimilation tool PESTPP-DA, which is part of the PEST++ framework (White et al., 2020), run under the pyEMU environment (White et al., 2016).

## 3. Application of remote sensing in integrated hydrological models

RS can contribute to IHMs by assigning: i) climate forcings, also
## Integrated hydrological modelling flowchart

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Flowchart summarizing the use of remote sensing and hydrogeophysics in integrated hydrological models.</figcaption></figure>

referred to as forcings (Section 3.1); ii) parameters of land surface and shallow soil (Section 3.2); iii) boundary conditions (Section 3.3); and iv) observations to constrain calibrated state variables of an IHM (Section 3.4). Climate forcings, system parameters and boundary conditions determine IHM solutions, while observations (RS, hydrogeophysics and in-situ) help to constrain simulated state variables in IHM calibration and data assimilation (Fig. 1).

### 3.1. Climate forcings

The common IHM climate forcings include effective precipitation and potential evapotranspiration (PET). Effective precipitation is the difference between precipitation and rainfall interception by plants, or shortly interception. Certain IHMs, which do not include inherent estimate of plant interception, e.g. MODFLOW 6, similarly require effective PET (PET minus interception) climate forcing. All the climate forcings' components are spatiotemporally variable, and as such, can be well-captured by RS, being particularly useful if bias-corrected by ground measurements.

#### 3.1.1. Precipitation

Precipitation is the most important climate forcing input of IHMs. The standard way of obtaining precipitation measurements is by tipping bucket rain gauges connected with logging devices. However, the representativeness of such measuring devices is restricted to an area of between 250 and 3,000 m² only (Kidd et al., 2017), depending on terrain complexity. Hence, considering the large spatiotemporal rainfall variability, particularly in topographically complex, water-limited environments (Derin and Yilmaz, 2014; Gebremedhin et al., 2021; Rahmawati and Lubczynski, 2018), an unrealistically large number of rain gauges would be required to assess rainfall properly, for example by interpolation (Camera et al., 2014). However, even that does not guarantee a reliable precipitation retrieval, as the high spatiotemporal variability of precipitation makes it rather impossible to capture precipitation reliably by interpolation. In this regard, earth observation of

precipitation provides a valuable contribution to spatiotemporal precipitation assessment.

Background information on the remote sensing of precipitation can be found in Levizzani and Cattani (2019) and Tang et al. (2020). Many RS precipitation products are freely available (Le Coz and van de Giesen, 2020) and do not require specialised RS knowledge. Sun et al., (2018a) compared 30 satellite precipitation products with global or reanalysis precipitation data sets. They found large differences in the magnitude and in the variability of the precipitation estimates and also that the spatial coverage of surface stations primarily constrained the reliability of precipitation estimates. Although the reliability (accuracy, spatial resolution, temporal resolution) of satellite precipitation products is continuously improved, they are all affected by biases of various types and origins. Therefore, the final reliability of precipitation estimates still strongly depends on: (i) the number and spatial distribution of gauges for bias correction, guidelines on the optimal location of rain gauges can be found in Morsy et al. (2021) and Di et al. (2020); (ii) the method of bias correction (see below); and also (iii) the accuracy and spatiotemporal resolution of a selected RS precipitation product.

There are many methods of rainfall bias correction (Ghimire et al., 2019). The accuracy of such methods, for a given number of gauges, their spatial distribution relative to the spatial scale of an IHM assessment and required by IHM temporal distribution, depends on how the bias is distributed (spatiotemporally). That distribution is influenced by various natural factors contributing differently to the bias (Gebremedhin et al., 2021; Rahmawati and Lubczynski, 2018). The most influencing bias correction factor is topography. Because rainfall increases with altitude, the bias increases, too. Hence, simple bias correction methods, such as those applied in the flat Kalahari Desert by Lekula et al. (2018), are not applicable in topographically complex areas where more sophisticated methods of merging satellite and daily gauge precipitation data are required; for example, the geographically weighted regression method is one of such methods, successfully applied for the bias correction of daily MPEG (Multi-sensor Precipitation Estimate-Geostationary), ~3 km spatial and 15 min temporal resolution and
## Integrated hydrological modelling flowchart

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Flowchart summarizing the use of remote sensing and hydrogeophysics in integrated hydrological models.</figcaption></figure>

referred to as forcings (Section 3.1); ii) parameters of land surface and shallow soil (Section 3.2); iii) boundary conditions (Section 3.3); and iv) observations to constrain calibrated state variables of an IHM (Section 3.4). Climate forcings, system parameters and boundary conditions determine IHM solutions, while observations (RS, hydrogeophysics and in-situ) help to constrain simulated state variables in IHM calibration and data assimilation (Fig. 1).

### 3.1. Climate forcings

The common IHM climate forcings include effective precipitation and potential evapotranspiration (PET). Effective precipitation is the difference between precipitation and rainfall interception by plants, or shortly interception. Certain IHMs, which do not include inherent estimate of plant interception, e.g. MODFLOW 6, similarly require effective PET (PET minus interception) climate forcing. All the climate forcings' components are spatiotemporally variable, and as such, can be well-captured by RS, being particularly useful if bias-corrected by ground measurements.

#### 3.1.1. Precipitation

Precipitation is the most important climate forcing input of IHMs. The standard way of obtaining precipitation measurements is by tipping bucket rain gauges connected with logging devices. However, the representativeness of such measuring devices is restricted to an area of between 250 and 3,000 m² only (Kidd et al., 2017), depending on terrain complexity. Hence, considering the large spatiotemporal rainfall variability, particularly in topographically complex, water-limited environments (Derin and Yilmaz, 2014; Gebremedhin et al., 2021; Rahmawati and Lubczynski, 2018), an unrealistically large number of rain gauges would be required to assess rainfall properly, for example by interpolation (Camera et al., 2014). However, even that does not guarantee a reliable precipitation retrieval, as the high spatiotemporal variability of precipitation makes it rather impossible to capture precipitation reliably by interpolation. In this regard, earth observation of

precipitation provides a valuable contribution to spatiotemporal precipitation assessment.

Background information on the remote sensing of precipitation can be found in Levizzani and Cattani (2019) and Tang et al. (2020). Many RS precipitation products are freely available (Le Coz and van de Giesen, 2020) and do not require specialised RS knowledge. Sun et al., (2018a) compared 30 satellite precipitation products with global or reanalysis precipitation data sets. They found large differences in the magnitude and in the variability of the precipitation estimates and also that the spatial coverage of surface stations primarily constrained the reliability of precipitation estimates. Although the reliability (accuracy, spatial resolution, temporal resolution) of satellite precipitation products is continuously improved, they are all affected by biases of various types and origins. Therefore, the final reliability of precipitation estimates still strongly depends on: (i) the number and spatial distribution of gauges for bias correction, guidelines on the optimal location of rain gauges can be found in Morsy et al. (2021) and Di et al. (2020); (ii) the method of bias correction (see below); and also (iii) the accuracy and spatiotemporal resolution of a selected RS precipitation product.

There are many methods of rainfall bias correction (Ghimire et al., 2019). The accuracy of such methods, for a given number of gauges, their spatial distribution relative to the spatial scale of an IHM assessment and required by IHM temporal distribution, depends on how the bias is distributed (spatiotemporally). That distribution is influenced by various natural factors contributing differently to the bias (Gebremedhin et al., 2021; Rahmawati and Lubczynski, 2018). The most influencing bias correction factor is topography. Because rainfall increases with altitude, the bias increases, too. Hence, simple bias correction methods, such as those applied in the flat Kalahari Desert by Lekula et al. (2018), are not applicable in topographically complex areas where more sophisticated methods of merging satellite and daily gauge precipitation data are required; for example, the geographically weighted regression method is one of such methods, successfully applied for the bias correction of daily MPEG (Multi-sensor Precipitation Estimate-Geostationary), ~3 km spatial and 15 min temporal resolution and
CHIRPS (Climate Hazards Group InfraRed Precipitation with Station), ~5 km spatial and 1-day temporal resolution satellite rainfall products in the topographically complex Upper Tekeze Basin in Ethiopia (Gebremedhin et al., 2021). Bias correction is also dependent on the scale of the temporal assessment. At coarse temporal scales, i.e. monthly, decadal or even weekly, the temporal data accumulation 'conveniently' reduces the systematic bias (AghaKouchak et al., 2012; Tian et al., 2009), but the daily bias correction is by far more demanding.

An optimal selection of a satellite precipitation product requires choosing the most accurately performing product at a given study area, considering its spatiotemporal characteristics and matching it with the IHM grid size and stress period. In the scientific literature, dozens of papers validate various satellite precipitation products at different locations to define the optimal one for an investigated area. However, most of the precipitation products used in those studies are at relatively coarse spatiotemporal resolution. To our knowledge, currently, the highest spatial resolution of satellite rainfall products, with at least daily temporal resolution (typically required by IHMs), are provided by MPEG and CHIRPS. The difference between the two, except for the spatiotemporal resolution differences (see above), is that CHIPRS, in contrast to MPEG, is inherently bias-corrected by selected gauging stations. Among coarser products, the IMERG (Integrated Multisatellite Retrievals for GPM) algorithm combining TRMM (Tropical Rainfall Measuring Mission) and GPM (Global Precipitation Measurement) with spatial resolution 0.1° (~10x10 km) and temporal 30 min, seems to be most interesting at present. IMERG appropriately detects and estimates regional precipitation patterns and reproduces extreme precipitation better than other coarse satellite products, although its performance in complex, especially mountainous terrains, needs improvement (Pradhan et al., 2022). Another global precipitation product MSWEP (Multi-Source Weighted-Ensemble Precipitation) merging gauge, satellite, and reanalysis data (Beck et al., 2019), is available since 1979 to real-time and has also a spatial resolution of 0.1° but a coarser temporal resolution of 3 h.

### 3.1.2. Rainfall interception loss by plants

Rainfall interception loss by plants, hereafter referred shortly as interception, represents part of rainfall never reaching the ground surface, so it is unavailable for infiltration or runoff. It is a major component of surface evaporation (Daoud et al., 2022) and is part of ET. In vegetated areas, interception can represent a considerable portion of ET. As interception reduces rainfall to effective precipitation and PET to effective PET, disregarding it or its erroneous estimation automatically introduces errors in the subsequent model-simulated processes (Save-nije, 2004). Experimentally, the interception of trees is measured by subtracting the sum of measured throughfall and stemflow from rainfall measurement (Ghimire et al., 2017; Ghimire et al., 2012), while cumbersome lysimeters (Crouse et al., 1966; McMillan and Burgy, 1960) or specific devices such as interception tubes (Demir et al., 2022), can estimate the interception of shorter plants such as grasses. However, in large areas with various land cover types, field experimental assessment of interception is not feasible, but IHMs still require its spatiotemporal estimate.

The most direct RS-based spatiotemporal interception assessment involves upscaling of field interception measurements by applying very high-resolution satellite images, as proposed by Hassan et al. (2017), who scaled field tree interception measurements by very fine (~1 m) spatial resolution IKONOS and WORLDVIEW II satellite images. Another option is to use analytical interception models (e.g. Gash, Rutter or Liu), evaluated by Eliades et al. (2022). These models were developed to quantify the temporal variability of interception of single trees based on field experimental measurements. Still, with the help of RS, they can be adapted to spatiotemporal assessment (Cui et al., 2015), because not only their hydrometeorological spatiotemporal variables, but also biophysical spatiotemporal variables (fraction of canopy cover and the canopy storage capacity) are definable from RS Leaf Area Index

estimates (Vegas Galdos et al., 2012; Zheng and Moskal, 2009). There are various RS adaptations of interception models, especially of the Gash model (Gash et al., 1995), e.g. at the global scale by Miralles et al. (2010) and at the catchment scale by Cui and Jia (2014); Gebremedhin et al. (2023). An example of spatiotemporal RS-based Gash model implementation in a catchment scale IHM (MODFLOW 6) is presented by Daoud et al. (2022).

Readily available interception products are still scarce. To our knowledge, the only such product is the WaPOR from the FAO portal at 250 m spatial resolution and decadal temporal resolution. However, WaPOR is available only for a selected part of Africa. Besides, it should be validated before it is used because compared to the RS Gash model, the WaPOR indicated substantial underestimation of interception and unrealistically low and invariable seasonal dynamics, likely because it simplifies the interception calculation, neglecting plant canopy storage capacity (Gebremedhin et al., 2023).

### 3.1.3. Potential evapotranspiration

PET is a controversial term in hydrology (Gebremedhin et al., 2022; Miralles et al., 2020; Savenije, 2004), as it is based on different approaches and definitions. However, regardless of the terminology, PET is a useful concept and an important IHM climate forcing, defined as 'the maximum amount of water that can be evaporated from soil, from surface features and from the subsurface as well as transpired by plants, under unlimited water availability condition' (Gebremedhin et al., 2022). The most common way of estimating PET (Allen et al., 1998) is as the product of the FAO reference evapotranspiration (ETo) and crop factor (Kc). As the FAO Kc refers exclusively to crops, Gebremedhin et al. (2022) addressed Kc as a land use-land cover factor, which refers not only to crops but also to natural vegetation.

FEWSNET (~111 km grid) and DMETREF (3 km grid) portals claim to provide RS-based, daily ETo corresponding with the FAO ETo. Considering the spatial resolution, the latter is more interesting, although it is affected by heat advection bias (Gebremedhin et al., 2022). Heat advection depends on the near-surface air temperature (Trigo et al., 2018). Hence, the bias can be corrected as a function of temperature. Gebremedhin et al. (2022) corrected ETo by applying ERA-5-Land temperature product (corresponding with measurements at the 2 m height), merging it with in-situ temperature measurements at the 2 m height, and applying the geographically weighted regression method.

For crops, the Kc is already well-defined by lysimetric measurements (Allen et al., 1998), but not for most natural land covers. The lysimetric measurements are cumbersome and inefficient, so they are generally not applicable for large area assessment, only for Kc calibration and validation. For spatiotemporal Kc variability in large areas, RS provides a valuable solution as Kc is linearly dependent on vegetation indices, particularly on NDVI (Rafn et al., 2008), and at given climatic conditions can be replicated from one area to another (Campos et al., 2010; Campos et al., 2013; Choudhury et al., 1994). An overview of recent improvements in applying RS to map Kc is presented by Pôças et al. (2020).

The Moderate Resolution Imaging Spectroradiometer (MODIS) provides an RS-based PET product at 500 m spatial and 8-day temporal resolution. However, the input data for this product is derived from the Global Modeling and Assimilation Office meteorological data at 1° x 1.25° (~111 km x 138.75 km) spatial resolution applying spatial smoothing (Mu et al., 2007). While this is acceptable for global or continental-scale IHMs, at the catchment-scale, it may lead to substantial biases (Mu et al., 2011). Another RS product, GLEAM, provides PET at daily resolution but at the global spatial resolution of 0.25° and is based on the Priestley and Taylor equation driven by observations of surface net radiation and near-surface air temperature (Martens et al., 2017).

### 3.1.4. Implementation of climate forcings in integrated hydrological models

IHM codes differ in their requirements for climate forcings as input. For example, MODFLOW 6 calculates evapotranspiration from the
CHIRPS (Climate Hazards Group InfraRed Precipitation with Station), ~5 km spatial and 1-day temporal resolution satellite rainfall products in the topographically complex Upper Tekeze Basin in Ethiopia (Gebremedhin et al., 2021). Bias correction is also dependent on the scale of the temporal assessment. At coarse temporal scales, i.e. monthly, decadal or even weekly, the temporal data accumulation 'conveniently' reduces the systematic bias (AghaKouchak et al., 2012; Tian et al., 2009), but the daily bias correction is by far more demanding.

An optimal selection of a satellite precipitation product requires choosing the most accurately performing product at a given study area, considering its spatiotemporal characteristics and matching it with the IHM grid size and stress period. In the scientific literature, dozens of papers validate various satellite precipitation products at different locations to define the optimal one for an investigated area. However, most of the precipitation products used in those studies are at relatively coarse spatiotemporal resolution. To our knowledge, currently, the highest spatial resolution of satellite rainfall products, with at least daily temporal resolution (typically required by IHMs), are provided by MPEG and CHIRPS. The difference between the two, except for the spatiotemporal resolution differences (see above), is that CHIPRS, in contrast to MPEG, is inherently bias-corrected by selected gauging stations. Among coarser products, the IMERG (Integrated Multisatellite Retrievals for GPM) algorithm combining TRMM (Tropical Rainfall Measuring Mission) and GPM (Global Precipitation Measurement) with spatial resolution 0.1° (~10x10 km) and temporal 30 min, seems to be most interesting at present. IMERG appropriately detects and estimates regional precipitation patterns and reproduces extreme precipitation better than other coarse satellite products, although its performance in complex, especially mountainous terrains, needs improvement (Pradhan et al., 2022). Another global precipitation product MSWEP (Multi-Source Weighted-Ensemble Precipitation) merging gauge, satellite, and reanalysis data (Beck et al., 2019), is available since 1979 to real-time and has also a spatial resolution of 0.1° but a coarser temporal resolution of 3 h.

### 3.1.2. Rainfall interception loss by plants

Rainfall interception loss by plants, hereafter referred shortly as interception, represents part of rainfall never reaching the ground surface, so it is unavailable for infiltration or runoff. It is a major component of surface evaporation (Daoud et al., 2022) and is part of ET. In vegetated areas, interception can represent a considerable portion of ET. As interception reduces rainfall to effective precipitation and PET to effective PET, disregarding it or its erroneous estimation automatically introduces errors in the subsequent model-simulated processes (Save-nije, 2004). Experimentally, the interception of trees is measured by subtracting the sum of measured throughfall and stemflow from rainfall measurement (Ghimire et al., 2017; Ghimire et al., 2012), while cumbersome lysimeters (Crouse et al., 1966; McMillan and Burgy, 1960) or specific devices such as interception tubes (Demir et al., 2022), can estimate the interception of shorter plants such as grasses. However, in large areas with various land cover types, field experimental assessment of interception is not feasible, but IHMs still require its spatiotemporal estimate.

The most direct RS-based spatiotemporal interception assessment involves upscaling of field interception measurements by applying very high-resolution satellite images, as proposed by Hassan et al. (2017), who scaled field tree interception measurements by very fine (~1 m) spatial resolution IKONOS and WORLDVIEW II satellite images. Another option is to use analytical interception models (e.g. Gash, Rutter or Liu), evaluated by Eliades et al. (2022). These models were developed to quantify the temporal variability of interception of single trees based on field experimental measurements. Still, with the help of RS, they can be adapted to spatiotemporal assessment (Cui et al., 2015), because not only their hydrometeorological spatiotemporal variables, but also biophysical spatiotemporal variables (fraction of canopy cover and the canopy storage capacity) are definable from RS Leaf Area Index

estimates (Vegas Galdos et al., 2012; Zheng and Moskal, 2009). There are various RS adaptations of interception models, especially of the Gash model (Gash et al., 1995), e.g. at the global scale by Miralles et al. (2010) and at the catchment scale by Cui and Jia (2014); Gebremedhin et al. (2023). An example of spatiotemporal RS-based Gash model implementation in a catchment scale IHM (MODFLOW 6) is presented by Daoud et al. (2022).

Readily available interception products are still scarce. To our knowledge, the only such product is the WaPOR from the FAO portal at 250 m spatial resolution and decadal temporal resolution. However, WaPOR is available only for a selected part of Africa. Besides, it should be validated before it is used because compared to the RS Gash model, the WaPOR indicated substantial underestimation of interception and unrealistically low and invariable seasonal dynamics, likely because it simplifies the interception calculation, neglecting plant canopy storage capacity (Gebremedhin et al., 2023).

### 3.1.3. Potential evapotranspiration

PET is a controversial term in hydrology (Gebremedhin et al., 2022; Miralles et al., 2020; Savenije, 2004), as it is based on different approaches and definitions. However, regardless of the terminology, PET is a useful concept and an important IHM climate forcing, defined as 'the maximum amount of water that can be evaporated from soil, from surface features and from the subsurface as well as transpired by plants, under unlimited water availability condition' (Gebremedhin et al., 2022). The most common way of estimating PET (Allen et al., 1998) is as the product of the FAO reference evapotranspiration (ETo) and crop factor (Kc). As the FAO Kc refers exclusively to crops, Gebremedhin et al. (2022) addressed Kc as a land use-land cover factor, which refers not only to crops but also to natural vegetation.

FEWSNET (~111 km grid) and DMETREF (3 km grid) portals claim to provide RS-based, daily ETo corresponding with the FAO ETo. Considering the spatial resolution, the latter is more interesting, although it is affected by heat advection bias (Gebremedhin et al., 2022). Heat advection depends on the near-surface air temperature (Trigo et al., 2018). Hence, the bias can be corrected as a function of temperature. Gebremedhin et al. (2022) corrected ETo by applying ERA-5-Land temperature product (corresponding with measurements at the 2 m height), merging it with in-situ temperature measurements at the 2 m height, and applying the geographically weighted regression method.

For crops, the Kc is already well-defined by lysimetric measurements (Allen et al., 1998), but not for most natural land covers. The lysimetric measurements are cumbersome and inefficient, so they are generally not applicable for large area assessment, only for Kc calibration and validation. For spatiotemporal Kc variability in large areas, RS provides a valuable solution as Kc is linearly dependent on vegetation indices, particularly on NDVI (Rafn et al., 2008), and at given climatic conditions can be replicated from one area to another (Campos et al., 2010; Campos et al., 2013; Choudhury et al., 1994). An overview of recent improvements in applying RS to map Kc is presented by Pôças et al. (2020).

The Moderate Resolution Imaging Spectroradiometer (MODIS) provides an RS-based PET product at 500 m spatial and 8-day temporal resolution. However, the input data for this product is derived from the Global Modeling and Assimilation Office meteorological data at 1° x 1.25° (~111 km x 138.75 km) spatial resolution applying spatial smoothing (Mu et al., 2007). While this is acceptable for global or continental-scale IHMs, at the catchment-scale, it may lead to substantial biases (Mu et al., 2011). Another RS product, GLEAM, provides PET at daily resolution but at the global spatial resolution of 0.25° and is based on the Priestley and Taylor equation driven by observations of surface net radiation and near-surface air temperature (Martens et al., 2017).

### 3.1.4. Implementation of climate forcings in integrated hydrological models

IHM codes differ in their requirements for climate forcings as input. For example, MODFLOW 6 calculates evapotranspiration from the
subsurface only, i.e. without considering interception. Consequently, when creating climate forcings' inputs, interception has to be externally subtracted from rainfall and from PET. In other IHM codes, where interception is internally calculated (e.g. GSFLOW), precipitation and PET can be directly used as climate forcing inputs.

If the time series of RS climate forcings contain data gaps (e.g. missing days in satellite products), then the missing data must be filled in, e.g. with statistical techniques. Moreover, RS climate forcings typically provide coarser spatial resolutions than the spatial resolution commonly used by IHMs. In such cases, the RS products can be down-scaled to the required IHM grid size. A recent review of downscaling methods of satellite-based precipitation estimates is provided by Abdollahipour et al. (2022). Conversely, if a satellite product has a finer resolution than the model grid, it can be aggregated, e.g. as a mean, to match with the IHM grid.

## 3.2. Parameters

In contrast to climate forcings and state variables, model parameters are time-invariant (at least at the temporal scale of IHM). Satellite RS has limited contributions to system parameterisation because RS signals can only inform land surface and very shallow subsurface, typically in the order of a few centimetres depth only.

Optical and multispectral (e.g. Landsat, Sentinel, etc.) satellite images provide georeferenced background for setting numerical model grids and for assignment of model features such as rivers, lakes, trees, urban areas, etc. (Qin and Liu, 2022; Tamiminia et al., 2020). They also provide landscape information, allowing land use-land cover classification.

A very important RS input for IHMs is a digital elevation model (DEM), which provides fundamental depictions of the three-dimensional shape of the Earth's surface (Guth et al., 2021). DEM provides important slope parameters influencing runoff and can help delineate watershed divides, river hydraulic gradient, etc. A DEM can also be used to study surface relief changes; for example, a combination of optical satellite images and the Shuttle Radar Topography Mission (SRTM) DEM was applied to map large-scale paleo-features, including a mega paleolake Chad, for the standalone groundwater model of the Lake Chad Basin (Leblanc et al., 2006a; Leblanc et al., 2007; Leblanc et al., 2006b).

RS DEMs are based on radar sensors (e.g. SRTM or ALOS), digital photogrammetry (IKONOS, ASTER) and Light Detection And Ranging (LiDAR) measurements. Hawker et al. (2019) list and compare global DEM products (coverage, acquisition years, sensor type, wavelength, resolution and vertical accuracy), focusing on performance comparison between SRTM, MERIT and TanDEM-X 90 applying airborne LiDAR as the reference. For IHMs, the most important element of a DEM is its vertical accuracy because it is used for referencing hydraulic heads and stages (for more information see Section 3.3). The vertical accuracy of different DEM products has been validated at various world locations; Shikoku Island (Japan) by Pakoksung and Takagi (2021); Maqu catchment of Tibetan Plateau by Li et al. (2021); and throughout the whole of Mexico by Carrera-Hernández (2021). Recently, a global 30 m resolution Copernicus surface DEM was presented for the first time, with buildings and forests removed by machine learning (Hawker et al., 2022).

One of the main RS applications for IHMs is RS-based soil parameterisation, i.e. mainly soil type and texture. General overviews of techniques, data sources and RS assessments of soil at different scales are available from Chen et al. (2022), Dewitte et al. (2012) and Mulder et al. (2011). Soil type and texture can be investigated by applying different RS data, such as: i) a combination of radar and optical data from Sentinel 1 and 2 (Bousbih et al., 2019); ii) multitemporal Landsat 8 data (Duan et al., 2022); or iii) hyperspectral data (Vibhute et al., 2015). Other IHM-important soil parameters can also be assessed through a combination of RS and other data types, mainly field-acquired data. For example, Honarbakhsh et al. (2022) used various RS data, including relief data in

combination with field-acquired, easy-to-measure soil properties, to determine spatially distributed soil-saturated hydraulic conductivity and pedo-transfer functions. Francés and Lubczynski (2011) used RS image processing of QuickBird, ASTER-GDEM (satellite images), aerial photographs, surface geophysics (EM-31), and reference soil samplings in geostatistical mixed linear model, all to determine the spatial variability of the clayey topsoil thickness for a catchment of ~ 19 km² at a spatial resolution of 25 m.

## 3.3. Boundary conditions

Boundary conditions are critical for model performance and should be defined with particular care. Adequate model boundaries align with physical boundaries (Anderson et al., 2015). Different types of boundary conditions can be defined on the basis of RS products. The lateral no-flow boundary conditions, such as watershed divides, can be determined from RS DEMs, while the contacts between different rock types, e.g. between permeable and impermeable rocks, as well as identification of fault lines or dyke locations, can be supported by digital processing of multispectral images. Optical and thermal remote sensing can also be used to detect and map groundwater discharge areas (Sass et al., 2014; Tweed et al., 2007).

Surface water bodies (e.g. lakes, rivers, etc.) are often assigned as boundary conditions of IHMs, especially if they are in hydraulic contact with groundwater. RS can identify surface water bodies and define their stages and related spatial extent. If they are assigned as specified head boundaries, then heads corresponding with these stages are forced, so all the uncertainty attributed to their assessment is carried over into the model results. The head-dependent boundary is a better solution as it is less constrained. In that solution, a temporally variable hydraulic gradient driving a flow (seepage) across the bottom of a surface water body is estimated from the difference between a stage observation and the simulated hydraulic head or simulated stage. For example, in MODFLOW, bottom seepages can be simulated through River or Reservoir Packages but also through more complex solutions offered by 'advanced' head-dependent boundaries, such as the SFR2 (for rivers) and Lake Packages of MODFLOW. In the latter two boundary cases, water exchanges between a surface water body and underlying unsaturated and saturated zone can be quantitatively controlled by observations of river discharges and lake stages, respectively. For example, El-Zehairy et al. (2018) used time series of river discharges and lake stage observations to constrain state variables in multivariate model calibration. More about RS observations of surface water body stages is presented in section 3.4.

Another type of boundary condition is groundwater abstraction. In irrigated areas, groundwater abstractions have been evaluated using RS of evapotranspiration (Ahmad et al., 2005; Vu et al., 2020). However, such RS assessment of groundwater withdrawals does not account for return flows to aquifers and does not differentiate between surface and groundwater, which can be problematic in the case of conjunctive water use for irrigation.

## 3.4. Remote sensing observations to constrain calibration of state variables

State variables represent in an IHM a simulated state of a hydrological system for a certain hydrological domain at a particular time and location. In calibration and data assimilation of IHMs, state variables are compared with observations (in-situ, hydrogeophysics or RS) to minimize the difference between them. Traditionally, groundwater, unsaturated zone and surface water flow model domains have been calibrated against classical in-situ observations, i.e. hydraulic heads, profile soil moisture (or matric potential) and surface water levels (stages) and discharges, respectively. In IHMs, various state variables of all model domains can be calibrated simultaneously. If the classical data are scarce, predictive uncertainty of IHMs can be reduced by calibration
subsurface only, i.e. without considering interception. Consequently, when creating climate forcings' inputs, interception has to be externally subtracted from rainfall and from PET. In other IHM codes, where interception is internally calculated (e.g. GSFLOW), precipitation and PET can be directly used as climate forcing inputs.

If the time series of RS climate forcings contain data gaps (e.g. missing days in satellite products), then the missing data must be filled in, e.g. with statistical techniques. Moreover, RS climate forcings typically provide coarser spatial resolutions than the spatial resolution commonly used by IHMs. In such cases, the RS products can be down-scaled to the required IHM grid size. A recent review of downscaling methods of satellite-based precipitation estimates is provided by Abdollahipour et al. (2022). Conversely, if a satellite product has a finer resolution than the model grid, it can be aggregated, e.g. as a mean, to match with the IHM grid.

## 3.2. Parameters

In contrast to climate forcings and state variables, model parameters are time-invariant (at least at the temporal scale of IHM). Satellite RS has limited contributions to system parameterisation because RS signals can only inform land surface and very shallow subsurface, typically in the order of a few centimetres depth only.

Optical and multispectral (e.g. Landsat, Sentinel, etc.) satellite images provide georeferenced background for setting numerical model grids and for assignment of model features such as rivers, lakes, trees, urban areas, etc. (Qin and Liu, 2022; Tamiminia et al., 2020). They also provide landscape information, allowing land use-land cover classification.

A very important RS input for IHMs is a digital elevation model (DEM), which provides fundamental depictions of the three-dimensional shape of the Earth's surface (Guth et al., 2021). DEM provides important slope parameters influencing runoff and can help delineate watershed divides, river hydraulic gradient, etc. A DEM can also be used to study surface relief changes; for example, a combination of optical satellite images and the Shuttle Radar Topography Mission (SRTM) DEM was applied to map large-scale paleo-features, including a mega paleolake Chad, for the standalone groundwater model of the Lake Chad Basin (Leblanc et al., 2006a; Leblanc et al., 2007; Leblanc et al., 2006b).

RS DEMs are based on radar sensors (e.g. SRTM or ALOS), digital photogrammetry (IKONOS, ASTER) and Light Detection And Ranging (LiDAR) measurements. Hawker et al. (2019) list and compare global DEM products (coverage, acquisition years, sensor type, wavelength, resolution and vertical accuracy), focusing on performance comparison between SRTM, MERIT and TanDEM-X 90 applying airborne LiDAR as the reference. For IHMs, the most important element of a DEM is its vertical accuracy because it is used for referencing hydraulic heads and stages (for more information see Section 3.3). The vertical accuracy of different DEM products has been validated at various world locations; Shikoku Island (Japan) by Pakoksung and Takagi (2021); Maqu catchment of Tibetan Plateau by Li et al. (2021); and throughout the whole of Mexico by Carrera-Hernández (2021). Recently, a global 30 m resolution Copernicus surface DEM was presented for the first time, with buildings and forests removed by machine learning (Hawker et al., 2022).

One of the main RS applications for IHMs is RS-based soil parameterisation, i.e. mainly soil type and texture. General overviews of techniques, data sources and RS assessments of soil at different scales are available from Chen et al. (2022), Dewitte et al. (2012) and Mulder et al. (2011). Soil type and texture can be investigated by applying different RS data, such as: i) a combination of radar and optical data from Sentinel 1 and 2 (Bousbih et al., 2019); ii) multitemporal Landsat 8 data (Duan et al., 2022); or iii) hyperspectral data (Vibhute et al., 2015). Other IHM-important soil parameters can also be assessed through a combination of RS and other data types, mainly field-acquired data. For example, Honarbakhsh et al. (2022) used various RS data, including relief data in

combination with field-acquired, easy-to-measure soil properties, to determine spatially distributed soil-saturated hydraulic conductivity and pedo-transfer functions. Francés and Lubczynski (2011) used RS image processing of QuickBird, ASTER-GDEM (satellite images), aerial photographs, surface geophysics (EM-31), and reference soil samplings in geostatistical mixed linear model, all to determine the spatial variability of the clayey topsoil thickness for a catchment of ~ 19 km² at a spatial resolution of 25 m.

## 3.3. Boundary conditions

Boundary conditions are critical for model performance and should be defined with particular care. Adequate model boundaries align with physical boundaries (Anderson et al., 2015). Different types of boundary conditions can be defined on the basis of RS products. The lateral no-flow boundary conditions, such as watershed divides, can be determined from RS DEMs, while the contacts between different rock types, e.g. between permeable and impermeable rocks, as well as identification of fault lines or dyke locations, can be supported by digital processing of multispectral images. Optical and thermal remote sensing can also be used to detect and map groundwater discharge areas (Sass et al., 2014; Tweed et al., 2007).

Surface water bodies (e.g. lakes, rivers, etc.) are often assigned as boundary conditions of IHMs, especially if they are in hydraulic contact with groundwater. RS can identify surface water bodies and define their stages and related spatial extent. If they are assigned as specified head boundaries, then heads corresponding with these stages are forced, so all the uncertainty attributed to their assessment is carried over into the model results. The head-dependent boundary is a better solution as it is less constrained. In that solution, a temporally variable hydraulic gradient driving a flow (seepage) across the bottom of a surface water body is estimated from the difference between a stage observation and the simulated hydraulic head or simulated stage. For example, in MODFLOW, bottom seepages can be simulated through River or Reservoir Packages but also through more complex solutions offered by 'advanced' head-dependent boundaries, such as the SFR2 (for rivers) and Lake Packages of MODFLOW. In the latter two boundary cases, water exchanges between a surface water body and underlying unsaturated and saturated zone can be quantitatively controlled by observations of river discharges and lake stages, respectively. For example, El-Zehairy et al. (2018) used time series of river discharges and lake stage observations to constrain state variables in multivariate model calibration. More about RS observations of surface water body stages is presented in section 3.4.

Another type of boundary condition is groundwater abstraction. In irrigated areas, groundwater abstractions have been evaluated using RS of evapotranspiration (Ahmad et al., 2005; Vu et al., 2020). However, such RS assessment of groundwater withdrawals does not account for return flows to aquifers and does not differentiate between surface and groundwater, which can be problematic in the case of conjunctive water use for irrigation.

## 3.4. Remote sensing observations to constrain calibration of state variables

State variables represent in an IHM a simulated state of a hydrological system for a certain hydrological domain at a particular time and location. In calibration and data assimilation of IHMs, state variables are compared with observations (in-situ, hydrogeophysics or RS) to minimize the difference between them. Traditionally, groundwater, unsaturated zone and surface water flow model domains have been calibrated against classical in-situ observations, i.e. hydraulic heads, profile soil moisture (or matric potential) and surface water levels (stages) and discharges, respectively. In IHMs, various state variables of all model domains can be calibrated simultaneously. If the classical data are scarce, predictive uncertainty of IHMs can be reduced by calibration
and/or data assimilation using unconventional observations (tracer concentrations, residence and travel times, exchange fluxes and temperature) as proposed by Schilling et al. (2019), but also by spatiotemporally variable RS (also hydrogeophysics) data originating from various satellite products.

### 3.4.1. Stages and discharges

RS-derived stages of surface water bodies, such as lakes, rivers etc., further referred shortly as stages, represent important observations for calibration of simulated stages or simulated stream discharges – in the latter case, the stream discharges are obtained from RS-derived stages. The easiest way to define the temporal variability of stages is by installing an automatic water level (stage) logging device. However, if such installation is not possible, RS provides the following possibilities to define stages of surface water bodies: i) indirect, by very high-resolution satellite observation of the extent of a surface water body, using an extent-stage relationship; ii) direct, by satellite altimetry. An example of the indirect RS estimate is how decadal time series of fluctuations for Lake Chad in Africa were inferred from high-resolution optical and thermal data (Leblanc et al., 2007; Leblanc et al., 2011). For smaller surface water bodies like rivers, the very high resolution (~1 m or better) optical RS monitoring showed a large potential for modelling ungauged catchments. For example, Sun et al., (2018b) used commercial, very high-resolution QuickBird, Ikonos and WorldView-1 images to calibrate a hydrological model of the upper Yalongjiang River Basin on the Tibetan Plateau. A similar study was proposed by Huang et al. (2022), who defined river stages by IKONOS, QuickBird, RapidEye-1/2/4/5, GeoEye-1, WorldView-2 and Landsat satellite images and applied them in a hydrological model to test its estimates of flows of the Lhasa, Salween, Mekong, and Yangtze Rivers at sections of ~100 m width. The reliability of the simulated discharges was confirmed by high Nash-Sutcliffe efficiency coefficients ranging between 0.5 and 0.8.

Direct RS measurement of stages of surface water bodies can be performed by high-resolution satellite altimetry, i.e. radar or laser altimetry. Satellite radar altimetry measures the time it takes for a radar pulse to travel from the satellite to the surface and back to the satellite receiver. Satellite laser altimetry does the same, only using light (photons) pulses. The Copernicus Sentinel 3A/3B, with a 300 m along-track and 1.64 km across-track footprint, was suitable for monitoring stages of wetlands in a direction parallel to the satellite tracks (Kittel et al., 2021). Deidda et al. (2021) showed, based on 22 Italian river sections with widths ranging from 50 to 555 m, that the stages of rivers can be quite accurately measured by free Sentinel 3A/3B and Jason 2/3 data. Garkoti and Kundapura (2021) used the same Copernicus altimetric missions to estimate river discharge solely using RS. The estimation was based on the simplified Manning's equation, assuming that a river's wetted perimeter was much larger than the river depth, which applies to discharge calculations of very wide rivers only. The river depth was estimated using satellite altimetry, the river bed slope by SRTM DEM and the river widths from the RivWidthCloud algorithm (Yang et al., 2020) processing of Sentinel 1 and Sentinel 2 images in Google Earth Engine. However, a cross-section geometry of a river is not always available, nor is the assumption of a wetted perimeter ≫ river depth, valid. For such cases, Jiang et al. (2021) proposed an approach based on a 1D MIKE HYDRO, i.e., a hydrodynamic model, combining cross-section geometry and roughness into one conveyance parameter. They used river widths from CryoSat-2 radar altimetry and Landsat imagery extracted using the RivWidthCloud algorithm in Google Earth Engine for model calibration.

For monitoring stages of large rivers wider than 100 m and for monitoring of surface water bodies, such as for example lakes, with a surface area larger than 250x250 m, very promising is the SWOT (Surface water and Ocean Topography) satellite hydrology mission, which was launched on 16 December 2022. This interferometry mission is expected to provide continuous information on the dynamics of surface

water bodies. One of its important products is SWORD (SWOT River Database), which contains data describing different river attributes such as network structure, water surface elevation, river width, and river type, so it can serve as a framework for global hydrological analyses using models, in situ measurements, and additional satellite observations (Altenau et al., 2021).

The disadvantage of RS radar altimetry in surface water body assessment is that it is still restricted to relatively large rivers, i.e. larger than ~100 m in width. The RS laser altimetry from ATLAS (Advanced Topographic Laser Altimeter System) onboard the ICESat-2 satellite mission, launched on 15 September 2018, overcomes this limitation. It offers a unique opportunity to measure river (Coppo Frias et al., 2023) and lake stages (Liu et al., 2022) at a much higher resolution than radar altimetry. For example, the along-track resolution of the ATL03 product is 0.7 m. The ICESat-2 not only allows for accurate assessment of water stages but also to map river (and lake) cross sections (Coppo Frias et al., 2023) and river water surface slopes (Scherer et al., 2022), both extremely useful for retrieval of river discharges even in non-uniform flow condition (Liu et al., 2023). A significant enhancement of surface water body monitoring is also expected through using ICESat-2 measurement via SWORD. For example, Christoffersen et al. (2023) developed an R package called ICE2WSS, which provides the opportunity to estimate the water surface slopes along the river centre lines in the SWORD using the ICESat-2 ATL13 water surface elevation measurements. They stated that the method works for slopes larger than 0.5 cm/km and river widths larger than 50 m. The main disadvantage of ICESat-2 is its limited data availability; as the mission was designed to monitor ice at the southern and northern poles, its data coverage declines towards the equator.

### 3.4.2. Evapotranspiration

RS estimates of ET represent useful observations capable of constraining simulated ET in model calibration and data assimilation. Originally, such applications used ET derived by RS solution of energy balance algorithms, such as the SEBAL (Bastiaanssen et al., 1998) applied to MODIS imagery to constrain the ET simulated by distributed rainfall-runoff SWAT (Arnold et al., 2012) model (not an IHM) of the Krishna basin (India) (Immerzeel and Droogers (2008). The study concluded that the RS-based ET could significantly improve the model calibration of the river discharge. A decade later, Wambura et al. (2018) also used RS-based ET as an observation to constrain a SWAT model, but instead of cumbersome and specialised energy balance ET calculations, they used a freely available MODIS ET product. In their study of the Wami River basin in Tanzania, a comparison was made between two cases: (a) using only a hydrograph as observation for model calibration; and (b) using hydrograph and MODIS-ET observations. They concluded that using the MODIS ET as an additional calibration constraint led to: (i) a better parameter calibration; (ii) a better capturing of the mean behaviour of ET; and (iii) a reduced prediction uncertainty.

The ET has also been used as an additional observation constraint in IHMs. For example, the MODIS ET product was used in the HydroGeo-Sphere IHM calibration next to groundwater levels and stream discharges of a hillslope catchment in Scotland by Ala-aho et al. (2017). The simulations revealed the critical role of groundwater in runoff generation but did not reproduce the remotely sensed evapotranspiration time series well. Hassan and Lubczynski (2024) also used MODIS ET product as additional (next to heads and 4 soil moisture profiles) constraint in GSFLOW calibration of a very fine, 5x5 m grid model, to investigate net recharge dependence on hydrotopes. Gelsinari et al. (2020) carried out a feasibility study of improving modelling results by applying a Kalman filter data assimilation of MODIS ET in a 1D conceptual unsaturated zone model coupled to MODFLOW 2005. The use of ET for data assimilation of the IHM was justified by the impact of the ET on groundwater levels and net recharge. The data assimilation improved not only the simulated evapotranspiration but also the net recharge, emphasizing the potential of RS ET as a calibration state variable.
and/or data assimilation using unconventional observations (tracer concentrations, residence and travel times, exchange fluxes and temperature) as proposed by Schilling et al. (2019), but also by spatiotemporally variable RS (also hydrogeophysics) data originating from various satellite products.

### 3.4.1. Stages and discharges

RS-derived stages of surface water bodies, such as lakes, rivers etc., further referred shortly as stages, represent important observations for calibration of simulated stages or simulated stream discharges – in the latter case, the stream discharges are obtained from RS-derived stages. The easiest way to define the temporal variability of stages is by installing an automatic water level (stage) logging device. However, if such installation is not possible, RS provides the following possibilities to define stages of surface water bodies: i) indirect, by very high-resolution satellite observation of the extent of a surface water body, using an extent-stage relationship; ii) direct, by satellite altimetry. An example of the indirect RS estimate is how decadal time series of fluctuations for Lake Chad in Africa were inferred from high-resolution optical and thermal data (Leblanc et al., 2007; Leblanc et al., 2011). For smaller surface water bodies like rivers, the very high resolution (~1 m or better) optical RS monitoring showed a large potential for modelling ungauged catchments. For example, Sun et al., (2018b) used commercial, very high-resolution QuickBird, Ikonos and WorldView-1 images to calibrate a hydrological model of the upper Yalongjiang River Basin on the Tibetan Plateau. A similar study was proposed by Huang et al. (2022), who defined river stages by IKONOS, QuickBird, RapidEye-1/2/4/5, GeoEye-1, WorldView-2 and Landsat satellite images and applied them in a hydrological model to test its estimates of flows of the Lhasa, Salween, Mekong, and Yangtze Rivers at sections of ~100 m width. The reliability of the simulated discharges was confirmed by high Nash-Sutcliffe efficiency coefficients ranging between 0.5 and 0.8.

Direct RS measurement of stages of surface water bodies can be performed by high-resolution satellite altimetry, i.e. radar or laser altimetry. Satellite radar altimetry measures the time it takes for a radar pulse to travel from the satellite to the surface and back to the satellite receiver. Satellite laser altimetry does the same, only using light (photons) pulses. The Copernicus Sentinel 3A/3B, with a 300 m along-track and 1.64 km across-track footprint, was suitable for monitoring stages of wetlands in a direction parallel to the satellite tracks (Kittel et al., 2021). Deidda et al. (2021) showed, based on 22 Italian river sections with widths ranging from 50 to 555 m, that the stages of rivers can be quite accurately measured by free Sentinel 3A/3B and Jason 2/3 data. Garkoti and Kundapura (2021) used the same Copernicus altimetric missions to estimate river discharge solely using RS. The estimation was based on the simplified Manning's equation, assuming that a river's wetted perimeter was much larger than the river depth, which applies to discharge calculations of very wide rivers only. The river depth was estimated using satellite altimetry, the river bed slope by SRTM DEM and the river widths from the RivWidthCloud algorithm (Yang et al., 2020) processing of Sentinel 1 and Sentinel 2 images in Google Earth Engine. However, a cross-section geometry of a river is not always available, nor is the assumption of a wetted perimeter ≫ river depth, valid. For such cases, Jiang et al. (2021) proposed an approach based on a 1D MIKE HYDRO, i.e., a hydrodynamic model, combining cross-section geometry and roughness into one conveyance parameter. They used river widths from CryoSat-2 radar altimetry and Landsat imagery extracted using the RivWidthCloud algorithm in Google Earth Engine for model calibration.

For monitoring stages of large rivers wider than 100 m and for monitoring of surface water bodies, such as for example lakes, with a surface area larger than 250x250 m, very promising is the SWOT (Surface water and Ocean Topography) satellite hydrology mission, which was launched on 16 December 2022. This interferometry mission is expected to provide continuous information on the dynamics of surface

water bodies. One of its important products is SWORD (SWOT River Database), which contains data describing different river attributes such as network structure, water surface elevation, river width, and river type, so it can serve as a framework for global hydrological analyses using models, in situ measurements, and additional satellite observations (Altenau et al., 2021).

The disadvantage of RS radar altimetry in surface water body assessment is that it is still restricted to relatively large rivers, i.e. larger than ~100 m in width. The RS laser altimetry from ATLAS (Advanced Topographic Laser Altimeter System) onboard the ICESat-2 satellite mission, launched on 15 September 2018, overcomes this limitation. It offers a unique opportunity to measure river (Coppo Frias et al., 2023) and lake stages (Liu et al., 2022) at a much higher resolution than radar altimetry. For example, the along-track resolution of the ATL03 product is 0.7 m. The ICESat-2 not only allows for accurate assessment of water stages but also to map river (and lake) cross sections (Coppo Frias et al., 2023) and river water surface slopes (Scherer et al., 2022), both extremely useful for retrieval of river discharges even in non-uniform flow condition (Liu et al., 2023). A significant enhancement of surface water body monitoring is also expected through using ICESat-2 measurement via SWORD. For example, Christoffersen et al. (2023) developed an R package called ICE2WSS, which provides the opportunity to estimate the water surface slopes along the river centre lines in the SWORD using the ICESat-2 ATL13 water surface elevation measurements. They stated that the method works for slopes larger than 0.5 cm/km and river widths larger than 50 m. The main disadvantage of ICESat-2 is its limited data availability; as the mission was designed to monitor ice at the southern and northern poles, its data coverage declines towards the equator.

### 3.4.2. Evapotranspiration

RS estimates of ET represent useful observations capable of constraining simulated ET in model calibration and data assimilation. Originally, such applications used ET derived by RS solution of energy balance algorithms, such as the SEBAL (Bastiaanssen et al., 1998) applied to MODIS imagery to constrain the ET simulated by distributed rainfall-runoff SWAT (Arnold et al., 2012) model (not an IHM) of the Krishna basin (India) (Immerzeel and Droogers (2008). The study concluded that the RS-based ET could significantly improve the model calibration of the river discharge. A decade later, Wambura et al. (2018) also used RS-based ET as an observation to constrain a SWAT model, but instead of cumbersome and specialised energy balance ET calculations, they used a freely available MODIS ET product. In their study of the Wami River basin in Tanzania, a comparison was made between two cases: (a) using only a hydrograph as observation for model calibration; and (b) using hydrograph and MODIS-ET observations. They concluded that using the MODIS ET as an additional calibration constraint led to: (i) a better parameter calibration; (ii) a better capturing of the mean behaviour of ET; and (iii) a reduced prediction uncertainty.

The ET has also been used as an additional observation constraint in IHMs. For example, the MODIS ET product was used in the HydroGeo-Sphere IHM calibration next to groundwater levels and stream discharges of a hillslope catchment in Scotland by Ala-aho et al. (2017). The simulations revealed the critical role of groundwater in runoff generation but did not reproduce the remotely sensed evapotranspiration time series well. Hassan and Lubczynski (2024) also used MODIS ET product as additional (next to heads and 4 soil moisture profiles) constraint in GSFLOW calibration of a very fine, 5x5 m grid model, to investigate net recharge dependence on hydrotopes. Gelsinari et al. (2020) carried out a feasibility study of improving modelling results by applying a Kalman filter data assimilation of MODIS ET in a 1D conceptual unsaturated zone model coupled to MODFLOW 2005. The use of ET for data assimilation of the IHM was justified by the impact of the ET on groundwater levels and net recharge. The data assimilation improved not only the simulated evapotranspiration but also the net recharge, emphasizing the potential of RS ET as a calibration state variable.
Recently, Gelsinari et al. (2021) compared also the two ET assimilation cases, i.e., the simple, 1D conceptual unsaturated zone model coupled to MODFLOW 2005 and the complex, Richards' equation 1D SWAP unsaturated zone model coupled to MODFLOW 2005. Remarkably, the data assimilations of the two coupled cases provided similar results.

### 3.4.3. Soil moisture

RS soil moisture (SM) has been used less frequently in IHMs as an observation to constrain model state variables than ET. Ridler et al. (2014) proposed one of the first studies assimilating satellite SMOS (Soil Moisture and Ocean Salinity) data in the MIKE SHE IHM. The RS SM observations substantially improved the simulated SM down to 25 cm depth but did not improve deeper SM at 50 cm depth, nor did the simulated river discharge improve. A more recent study by Zhao et al. (2021) assimilated SM of the SMAP (Soil Moisture Active Passive) product, first in the surface community land model (CLM) and then in the IHM coupling of the same CLM with the ParFlow (subsurface). In contrast to CLM alone, which considered only vertical flow, the CLM-ParFlow IHM considered both vertical and lateral flows, showing less bias and similar random errors but a slightly smaller correlation with the measurements. Nevertheless, the CLM-ParFlow reproduced better soil moisture dynamics, emphasizing the advantage of multi-domain IHMs compared to single-domain standalone models (Zhao et al., 2021). The highest spatial SM resolution (1 km) is provided by SENTINEL-1 (C-SAR) product (https://land.copernicus.eu/global/products/ssm) available daily in Europe since October 2014. For example, the applicability of SENTINEL-1 topsoil SM for hydrological simulation with HEC-HMS model has been tested in Keramianos basin in Crete by (Alexakis et al., 2017).

### 3.4.4. Change of terrestrial water storage

Another observation type that can be applied as a constraint of state variables in IHM calibration and data assimilation is the temporal change of terrestrial water storage (TWS) recorded by the Gravity Recovery and Climate Experiment (GRACE) mission since 2002 (Humphrey et al., 2023). However, isolating the groundwater component from the TWS remains a challenge. This is generally achieved by removing from TWS soil moisture, snow and surface water components, which are estimated, depending on the complexity of a basin, using a combination of in situ data, satellite data and models as proposed by Frappart and Ramillien (2018). Different solution was proposed by Andrew et al. (2017), who showed that wavelet decomposition might also be used to split total water storage into different components, including soil water and groundwater.

The validations of GRACE groundwater storage changes have been carried out by directly comparing them with groundwater level changes in wells and with models. For example, Rateb et al. (2020) compared groundwater storage changes from GRACE over 15 years (2002–2017) in 14 major U.S. aquifers, with groundwater-level monitoring data in ~23,000 wells and with regional and global hydrological and land surface models. For validation of GRACE groundwater storage changes, IHMs have also been used. For example, Seyoum and Milewski (2016) developed a MIKE SHE model of the Northern High Plains in the USA that was calibrated against independent, in-situ (heads, stream flows, soil moisture) observations. They found a good agreement between GRACE TWS and the TWS of the MIKE SHE IHM (simulating the entire terrestrial water cycle).

The usefulness of GRACE observations to control IHM state variables is unquestionable, as GRACE data can help to further constrain IHM parameters compared to in situ data alone. Sun et al. (2012) used GRACE groundwater storage change to calibrate a MODFLOW-based groundwater model of a regional aquifer in the USA, adjusting hydraulic conductivity, specific yield, and a diffuse recharge multiplier simultaneously. Hu and Jiao (2015) used GRACE groundwater storage change to calibrate the hydraulic conductivity of a large-scale FEFLOW model of the Qaidam Basin, China. They converted GRACE groundwater

storage into an average groundwater head for each model cell. Li et al. (2019) assimilated a state-of-the-art TWS GRACE product into NASA's Catchment Land Surface Model at the global scale, with the goal of generating groundwater storage time series that were useful for drought monitoring and other applications. Evaluation using in situ data from nearly 4,000 wells showed that GRACE data assimilation improved the simulation of groundwater at regional and point scales, with estimation errors reduced by 36% and 10% and correlation improved by 16% and 22%, respectively.

The main limitation of GRACE data is its very coarse spatial resolution (not better than 100x100 km), restricting its use to large basins or continental scales. However, interestingly, land surface and hydrological models can be used for downscaling GRACE, e.g. Vishwakarma et al. (2021). Other downscaling approaches, based on artificial intelligence, were also able to provide useful data for IHMs. For example, Jyolsna et al. (2021) established GRACE TWS correlation with various land surface and hydroclimatic variables in two commonly used machine learning algorithms (multi-linear regression and random forest models) to downscale groundwater storage changes.

### 3.4.5. Multi-mission observations

The combined use of multi-mission satellite products provides IHM input of different types and different spatiotemporal resolutions. The multi-mission satellite products can be applied as observations next to in-situ data for multivariate calibration or data assimilation. The multi-mission data better constrain model solutions than single-mission data, so it also better reduces model uncertainty. The use of multi-mission satellite products in a SWAT model of the Okavango Basin was proposed by Milzow et al. (2011), who used a combination of: i) satellite altimetry from the RLA (River and Lake Altimetry) product of Montfort University based on various altimetry missions; ii) soil moisture from ASAR (Advanced Synthetic Aperture Radar) and iii) total water storage from the GRACE satellite mission. The satellite-based and ground-based observations were calibrated through an automatic, multivariate procedure, which confirmed the benefit of using different, multi-mission observations in the model calibration, but also emphasised a loss of model accuracy in simulating catchment outflow due to errors in the RS-obtained precipitation forcing. Dembele et al. (2020) used multi-mission satellite products in a multivariate calibration framework of a conceptual, distributed mesoscale Hydrologic Model (mHM) (not IHM) for the gauged Volta River basin in West Africa. In the mHM calibration, they simultaneously incorporated streamflow and three satellite products, i. e., ET from GLEAM (Global Land Evaporation Amsterdam Model; https://www.gleam.eu), SM from ESA-CCI (European Space Agency Climate Change Initiative; https://www.climate.esa.int) and TWS from GRACE (www2.jpl.nasa.gov/grace/). The MODIS land surface temperature data set was used for model evaluation. The study confirmed and demonstrated benefits of using multiple satellite data sets in a model calibration.

Very few studies, especially those using physically based IHMs, simulate both the surface water and groundwater domains while integrating multi-mission satellite products as observations to constrain model calibration. Two different satellite-based products, i.e. ET of GLEAM and SM of ESA-CCI, were used by Lopez et al. (2017) in the calibration of a large scale (~38,000 km² and a 10 by 10 km grid) hydrological PCRaster GLOBal Water Balance (PCR-GLOBWB) model in Morocco. The PCR-GLOBWB is not a physically based IHM but a leaky-bucket type water balance model, which involves all model domains, including the groundwater system simulated as a combination of two conceptual reservoirs. The multivariate model calibration with ET and SM controlled by in-situ river discharge measurements, provided better results than using ET or SM separately, substantially reducing the equifinality problem. Gaur et al. (2022) used two RS-based observation products, the ET of MODIS and the SM of ESA-CCI, in a physically based MIKE SHE IHM of the 19,276 km² and 5 by 5 km grid over Subarnarekha catchment in Eastern India. The multivariate calibration focused on the
Recently, Gelsinari et al. (2021) compared also the two ET assimilation cases, i.e., the simple, 1D conceptual unsaturated zone model coupled to MODFLOW 2005 and the complex, Richards' equation 1D SWAP unsaturated zone model coupled to MODFLOW 2005. Remarkably, the data assimilations of the two coupled cases provided similar results.

### 3.4.3. Soil moisture

RS soil moisture (SM) has been used less frequently in IHMs as an observation to constrain model state variables than ET. Ridler et al. (2014) proposed one of the first studies assimilating satellite SMOS (Soil Moisture and Ocean Salinity) data in the MIKE SHE IHM. The RS SM observations substantially improved the simulated SM down to 25 cm depth but did not improve deeper SM at 50 cm depth, nor did the simulated river discharge improve. A more recent study by Zhao et al. (2021) assimilated SM of the SMAP (Soil Moisture Active Passive) product, first in the surface community land model (CLM) and then in the IHM coupling of the same CLM with the ParFlow (subsurface). In contrast to CLM alone, which considered only vertical flow, the CLM-ParFlow IHM considered both vertical and lateral flows, showing less bias and similar random errors but a slightly smaller correlation with the measurements. Nevertheless, the CLM-ParFlow reproduced better soil moisture dynamics, emphasizing the advantage of multi-domain IHMs compared to single-domain standalone models (Zhao et al., 2021). The highest spatial SM resolution (1 km) is provided by SENTINEL-1 (C-SAR) product (https://land.copernicus.eu/global/products/ssm) available daily in Europe since October 2014. For example, the applicability of SENTINEL-1 topsoil SM for hydrological simulation with HEC-HMS model has been tested in Keramianos basin in Crete by (Alexakis et al., 2017).

### 3.4.4. Change of terrestrial water storage

Another observation type that can be applied as a constraint of state variables in IHM calibration and data assimilation is the temporal change of terrestrial water storage (TWS) recorded by the Gravity Recovery and Climate Experiment (GRACE) mission since 2002 (Humphrey et al., 2023). However, isolating the groundwater component from the TWS remains a challenge. This is generally achieved by removing from TWS soil moisture, snow and surface water components, which are estimated, depending on the complexity of a basin, using a combination of in situ data, satellite data and models as proposed by Frappart and Ramillien (2018). Different solution was proposed by Andrew et al. (2017), who showed that wavelet decomposition might also be used to split total water storage into different components, including soil water and groundwater.

The validations of GRACE groundwater storage changes have been carried out by directly comparing them with groundwater level changes in wells and with models. For example, Rateb et al. (2020) compared groundwater storage changes from GRACE over 15 years (2002–2017) in 14 major U.S. aquifers, with groundwater-level monitoring data in ~23,000 wells and with regional and global hydrological and land surface models. For validation of GRACE groundwater storage changes, IHMs have also been used. For example, Seyoum and Milewski (2016) developed a MIKE SHE model of the Northern High Plains in the USA that was calibrated against independent, in-situ (heads, stream flows, soil moisture) observations. They found a good agreement between GRACE TWS and the TWS of the MIKE SHE IHM (simulating the entire terrestrial water cycle).

The usefulness of GRACE observations to control IHM state variables is unquestionable, as GRACE data can help to further constrain IHM parameters compared to in situ data alone. Sun et al. (2012) used GRACE groundwater storage change to calibrate a MODFLOW-based groundwater model of a regional aquifer in the USA, adjusting hydraulic conductivity, specific yield, and a diffuse recharge multiplier simultaneously. Hu and Jiao (2015) used GRACE groundwater storage change to calibrate the hydraulic conductivity of a large-scale FEFLOW model of the Qaidam Basin, China. They converted GRACE groundwater

storage into an average groundwater head for each model cell. Li et al. (2019) assimilated a state-of-the-art TWS GRACE product into NASA's Catchment Land Surface Model at the global scale, with the goal of generating groundwater storage time series that were useful for drought monitoring and other applications. Evaluation using in situ data from nearly 4,000 wells showed that GRACE data assimilation improved the simulation of groundwater at regional and point scales, with estimation errors reduced by 36% and 10% and correlation improved by 16% and 22%, respectively.

The main limitation of GRACE data is its very coarse spatial resolution (not better than 100x100 km), restricting its use to large basins or continental scales. However, interestingly, land surface and hydrological models can be used for downscaling GRACE, e.g. Vishwakarma et al. (2021). Other downscaling approaches, based on artificial intelligence, were also able to provide useful data for IHMs. For example, Jyolsna et al. (2021) established GRACE TWS correlation with various land surface and hydroclimatic variables in two commonly used machine learning algorithms (multi-linear regression and random forest models) to downscale groundwater storage changes.

### 3.4.5. Multi-mission observations

The combined use of multi-mission satellite products provides IHM input of different types and different spatiotemporal resolutions. The multi-mission satellite products can be applied as observations next to in-situ data for multivariate calibration or data assimilation. The multi-mission data better constrain model solutions than single-mission data, so it also better reduces model uncertainty. The use of multi-mission satellite products in a SWAT model of the Okavango Basin was proposed by Milzow et al. (2011), who used a combination of: i) satellite altimetry from the RLA (River and Lake Altimetry) product of Montfort University based on various altimetry missions; ii) soil moisture from ASAR (Advanced Synthetic Aperture Radar) and iii) total water storage from the GRACE satellite mission. The satellite-based and ground-based observations were calibrated through an automatic, multivariate procedure, which confirmed the benefit of using different, multi-mission observations in the model calibration, but also emphasised a loss of model accuracy in simulating catchment outflow due to errors in the RS-obtained precipitation forcing. Dembele et al. (2020) used multi-mission satellite products in a multivariate calibration framework of a conceptual, distributed mesoscale Hydrologic Model (mHM) (not IHM) for the gauged Volta River basin in West Africa. In the mHM calibration, they simultaneously incorporated streamflow and three satellite products, i. e., ET from GLEAM (Global Land Evaporation Amsterdam Model; https://www.gleam.eu), SM from ESA-CCI (European Space Agency Climate Change Initiative; https://www.climate.esa.int) and TWS from GRACE (www2.jpl.nasa.gov/grace/). The MODIS land surface temperature data set was used for model evaluation. The study confirmed and demonstrated benefits of using multiple satellite data sets in a model calibration.

Very few studies, especially those using physically based IHMs, simulate both the surface water and groundwater domains while integrating multi-mission satellite products as observations to constrain model calibration. Two different satellite-based products, i.e. ET of GLEAM and SM of ESA-CCI, were used by Lopez et al. (2017) in the calibration of a large scale (~38,000 km² and a 10 by 10 km grid) hydrological PCRaster GLOBal Water Balance (PCR-GLOBWB) model in Morocco. The PCR-GLOBWB is not a physically based IHM but a leaky-bucket type water balance model, which involves all model domains, including the groundwater system simulated as a combination of two conceptual reservoirs. The multivariate model calibration with ET and SM controlled by in-situ river discharge measurements, provided better results than using ET or SM separately, substantially reducing the equifinality problem. Gaur et al. (2022) used two RS-based observation products, the ET of MODIS and the SM of ESA-CCI, in a physically based MIKE SHE IHM of the 19,276 km² and 5 by 5 km grid over Subarnarekha catchment in Eastern India. The multivariate calibration focused on the
intercomparison of RS-based ET and SM patterns with the simulated patterns, applying three spatial performance metrics as diagnostic tools: i) joint empirical orthogonal functions; ii) fractional skill score; and iii) spatial efficiency. Finally, the uncertainty associated with these patterns, was assessed through the quantile regression technique, which confirmed that the ET and SM improved the water balance of the catchment. An example of the evaluation (not calibration) of the fully-coupled ParFlow-CLM (PF-CONUS v.1) IHM, developed by Maxwell et al. (2015), applying several RS products and thousands of in-situ data, is presented by O'Neill et al. (2021). The transient model was run hourly for four years at the continental scale of the USA, using a 1 by 1 km grid. The study showed the benefit of using remote sensing data for model evaluation.

The multivariate data assimilation of multi-mission satellite products was proposed by Khaki et al. (2020), who tested the global water balance W3RA model over two different large basins, the Murray-Darling and Mississippi. The RS products used as observations included: soil moisture from SMOS and AMSR-E (Advanced Microwave Scanning Radiometer-Earth Observing System), Leaf Area Index from AVHRR (Advanced Very-High Resolution Radiometer) and terrestrial water storage from GRACE. Although the W3RA model was not an IHM, the study showed that simultaneous assimilation of observations from multiple satellite products, combined with parameter estimation, strongly improved model predictive capability compared with single satellite products or state estimation alone. To our knowledge, no published study has applied multi-mission satellite products in multivariate data assimilation of physically based IHM, at least not by the time of the release of this review.

# 4. Application of hydrogeophysics in integrated hydrological models

The usefulness of geophysics, so also hydrogeophysics, in obtaining input data for IHMs, depends on the presence of contrast in physical properties of the subsurface that can be inferred from hydrogeophysical surveys and their analysis. The subsurface physical properties or quantities that can be measured during such surveys include: electrical properties (resistivity, permittivity, induced and spontaneous polarization), seismic velocities, radioactivity, gravity, earth magnetic field, or hydrogen proton precession in Nuclear Magnetic Resonance, etc.

Three types of hydrogeophysical surveys as dependent on data acquisition platforms, are discussed considering their potential input for IHMs: i) surface hydrogeophysics, i.e. surveys carried out from the ground surface ii) drone-borne hydrogeophysics, i.e. surveys carried out from remotely controlled drones; iii) airborne hydrogeophysics, i.e. surveys carried out from piloted aircraft, mainly airplanes and helicopters.

## 4.1. Surface hydrogeophysics

Six different surface hydrogeophysical methods are discussed: electrical resistivity (ER), electromagnetic (EM), ground penetrating radar (GPR), seismic reflection, terrestrial gravimetry and magnetic resonance sounding (MRS).

### 4.1.1. Electrical resistivity

ER techniques (tomography and profiling) focus on mapping the subsurface's electrical resistivity by injecting electrical current into the ground and measuring the resulting electrical potential in a set of electrodes. In resistive and moderately resistive environments, ER technique can reach large depths by increasing the electrodes' distance (large array size), although losing resolution with depth. The ER technique is less sensitive to environmental noise than other hydrogeophysical methods, but it has moderate operational efficiency because of time-consuming field installation and data acquisition (Francés et al., 2014). Besides, the ER method requires good electrical contact with the

subsurface to inject the current. Therefore, its use in dry soils and outcropping rocks can be problematic. Nevertheless, the ER techniques have a long-proven record of contributing relevant hydro(geo)logical information (Hübner et al., 2015). The equipment's reasonable cost and the measured parameter's dependence upon soil and rock composition, porosity, interstitial water salinity, saturation, etc., make the technique still attractive. Since it lacks selectivity for water, optimal exploitation assumes at least general local geology and hydrogeology knowledge.

### 4.1.2. Electromagnetic

EM surveys measure apparent electrical conductivity either in the frequency domain electromagnetic (FDEM) (Siemon, 2009) or in the time domain electromagnetic (TDEM), the latter also known as transient electromagnetics (TEM) (Christiansen et al., 2006). FDEM measures the signal amplitude and phase of an EM-induced field, while TDEM measures signal amplitude and decay time constant following an EM pulse induced by a transmitter. Both EM methods assume that some of the hydrostratigraphic units are at least moderately conductive. The field FDEM survey is very efficient as no galvanic contact with the ground is needed to assess the electrical properties of the subsurface and, therefore, is typically adapted to provide lateral changes in the conductivity of the shallow subsurface. Its depth of investigation depends on the range of frequencies, coil spacing, signal-to-noise (S/N) ratio and model used for data inversion (Francés and Lubczynski, 2011). For example, with a fixed 5 m coil spacing on a rigid frame enhancing S/N ratio, and using the Slingram method with the MAX-MIN 18S equipment with 8 frequencies from 444 Hz to 56 kHz, Frances et al. (2014) were able to reach 50 m depth penetration in weathered and fractured granites, while realizing a survey transect of 1620 m with data acquired every 60 m within ~ 5 h. The field TDEM survey, in its classical implementation, is less efficient than FDEM, as it involves on the ground setup of electrical cables of a transmitter and receiver coils, with the size proportional to the required depth of investigation. However, when used with large coil separations and frequencies adjusted to the local mineral composition of targets (deep aquifers, deep basement boundary, freshwater-saltwater interface, etc.), the TDEM can map the subsurface even down to several hundred meters (Frances et al., 2015). Recent advances in the TDEM resulted also in efficient, towed TDEM instrument implementation, capable of high-resolution 3D imaging of the subsurface down to a depth of 70 m (Auken et al., 2019).

### 4.1.3. Ground penetrating radar

GPR (also known as Georadar) instruments are equipped with transmitting and receiving antennas with varying frequencies, typically from 40 MHz to 2 GHz. The transmitting antenna sends electromagnetic waves, which, if reflected at interfaces (e.g., at boundaries of subsurface layers or bodies differing by dielectric permittivity), are then registered by the receiving antenna. The time-wise operational efficiency of GPR is even larger than that of the FDEM, as the GPR can work in continuous mode being towed by a vehicle (Mahmoudzadeh et al., 2012). Theoretically, in resistive environments, GPR has a large penetration depth (up to 30 m), but in conductive environments, its signal can be attenuated even at a very shallow subsurface, i.e. ~ 1 m. GPR has many different applications, but in hydrology, its main application is the assessment of soil moisture (Lambot et al., 2008; Tran et al., 2015; Wu et al., 2022). Also, the bathymetry of water bodies and characteristics of bottom sediments, such as sediment porosity, can be done by GPR (Sambuelli and Bava, 2012). However, the GPR bathymetry survey works well only in freshwater environments, but with increased electric conductivity and dielectric permittivity of water, its penetration depth decreases. An assessment of water table depth is also considered a very important GPR application. However, the anticipated water table response is not observed in the GPR data at many sites (Annan et al., 1991). For example, in favourable, highly resistive weathered granite survey conditions, while working at 200 MHz GPR frequency, Mahmoudzadeh et al. (2012) could detect water table depth only down to ~
intercomparison of RS-based ET and SM patterns with the simulated patterns, applying three spatial performance metrics as diagnostic tools: i) joint empirical orthogonal functions; ii) fractional skill score; and iii) spatial efficiency. Finally, the uncertainty associated with these patterns, was assessed through the quantile regression technique, which confirmed that the ET and SM improved the water balance of the catchment. An example of the evaluation (not calibration) of the fully-coupled ParFlow-CLM (PF-CONUS v.1) IHM, developed by Maxwell et al. (2015), applying several RS products and thousands of in-situ data, is presented by O'Neill et al. (2021). The transient model was run hourly for four years at the continental scale of the USA, using a 1 by 1 km grid. The study showed the benefit of using remote sensing data for model evaluation.

The multivariate data assimilation of multi-mission satellite products was proposed by Khaki et al. (2020), who tested the global water balance W3RA model over two different large basins, the Murray-Darling and Mississippi. The RS products used as observations included: soil moisture from SMOS and AMSR-E (Advanced Microwave Scanning Radiometer-Earth Observing System), Leaf Area Index from AVHRR (Advanced Very-High Resolution Radiometer) and terrestrial water storage from GRACE. Although the W3RA model was not an IHM, the study showed that simultaneous assimilation of observations from multiple satellite products, combined with parameter estimation, strongly improved model predictive capability compared with single satellite products or state estimation alone. To our knowledge, no published study has applied multi-mission satellite products in multivariate data assimilation of physically based IHM, at least not by the time of the release of this review.

# 4. Application of hydrogeophysics in integrated hydrological models

The usefulness of geophysics, so also hydrogeophysics, in obtaining input data for IHMs, depends on the presence of contrast in physical properties of the subsurface that can be inferred from hydrogeophysical surveys and their analysis. The subsurface physical properties or quantities that can be measured during such surveys include: electrical properties (resistivity, permittivity, induced and spontaneous polarization), seismic velocities, radioactivity, gravity, earth magnetic field, or hydrogen proton precession in Nuclear Magnetic Resonance, etc.

Three types of hydrogeophysical surveys as dependent on data acquisition platforms, are discussed considering their potential input for IHMs: i) surface hydrogeophysics, i.e. surveys carried out from the ground surface ii) drone-borne hydrogeophysics, i.e. surveys carried out from remotely controlled drones; iii) airborne hydrogeophysics, i.e. surveys carried out from piloted aircraft, mainly airplanes and helicopters.

## 4.1. Surface hydrogeophysics

Six different surface hydrogeophysical methods are discussed: electrical resistivity (ER), electromagnetic (EM), ground penetrating radar (GPR), seismic reflection, terrestrial gravimetry and magnetic resonance sounding (MRS).

### 4.1.1. Electrical resistivity

ER techniques (tomography and profiling) focus on mapping the subsurface's electrical resistivity by injecting electrical current into the ground and measuring the resulting electrical potential in a set of electrodes. In resistive and moderately resistive environments, ER technique can reach large depths by increasing the electrodes' distance (large array size), although losing resolution with depth. The ER technique is less sensitive to environmental noise than other hydrogeophysical methods, but it has moderate operational efficiency because of time-consuming field installation and data acquisition (Francés et al., 2014). Besides, the ER method requires good electrical contact with the

subsurface to inject the current. Therefore, its use in dry soils and outcropping rocks can be problematic. Nevertheless, the ER techniques have a long-proven record of contributing relevant hydro(geo)logical information (Hübner et al., 2015). The equipment's reasonable cost and the measured parameter's dependence upon soil and rock composition, porosity, interstitial water salinity, saturation, etc., make the technique still attractive. Since it lacks selectivity for water, optimal exploitation assumes at least general local geology and hydrogeology knowledge.

### 4.1.2. Electromagnetic

EM surveys measure apparent electrical conductivity either in the frequency domain electromagnetic (FDEM) (Siemon, 2009) or in the time domain electromagnetic (TDEM), the latter also known as transient electromagnetics (TEM) (Christiansen et al., 2006). FDEM measures the signal amplitude and phase of an EM-induced field, while TDEM measures signal amplitude and decay time constant following an EM pulse induced by a transmitter. Both EM methods assume that some of the hydrostratigraphic units are at least moderately conductive. The field FDEM survey is very efficient as no galvanic contact with the ground is needed to assess the electrical properties of the subsurface and, therefore, is typically adapted to provide lateral changes in the conductivity of the shallow subsurface. Its depth of investigation depends on the range of frequencies, coil spacing, signal-to-noise (S/N) ratio and model used for data inversion (Francés and Lubczynski, 2011). For example, with a fixed 5 m coil spacing on a rigid frame enhancing S/N ratio, and using the Slingram method with the MAX-MIN 18S equipment with 8 frequencies from 444 Hz to 56 kHz, Frances et al. (2014) were able to reach 50 m depth penetration in weathered and fractured granites, while realizing a survey transect of 1620 m with data acquired every 60 m within ~ 5 h. The field TDEM survey, in its classical implementation, is less efficient than FDEM, as it involves on the ground setup of electrical cables of a transmitter and receiver coils, with the size proportional to the required depth of investigation. However, when used with large coil separations and frequencies adjusted to the local mineral composition of targets (deep aquifers, deep basement boundary, freshwater-saltwater interface, etc.), the TDEM can map the subsurface even down to several hundred meters (Frances et al., 2015). Recent advances in the TDEM resulted also in efficient, towed TDEM instrument implementation, capable of high-resolution 3D imaging of the subsurface down to a depth of 70 m (Auken et al., 2019).

### 4.1.3. Ground penetrating radar

GPR (also known as Georadar) instruments are equipped with transmitting and receiving antennas with varying frequencies, typically from 40 MHz to 2 GHz. The transmitting antenna sends electromagnetic waves, which, if reflected at interfaces (e.g., at boundaries of subsurface layers or bodies differing by dielectric permittivity), are then registered by the receiving antenna. The time-wise operational efficiency of GPR is even larger than that of the FDEM, as the GPR can work in continuous mode being towed by a vehicle (Mahmoudzadeh et al., 2012). Theoretically, in resistive environments, GPR has a large penetration depth (up to 30 m), but in conductive environments, its signal can be attenuated even at a very shallow subsurface, i.e. ~ 1 m. GPR has many different applications, but in hydrology, its main application is the assessment of soil moisture (Lambot et al., 2008; Tran et al., 2015; Wu et al., 2022). Also, the bathymetry of water bodies and characteristics of bottom sediments, such as sediment porosity, can be done by GPR (Sambuelli and Bava, 2012). However, the GPR bathymetry survey works well only in freshwater environments, but with increased electric conductivity and dielectric permittivity of water, its penetration depth decreases. An assessment of water table depth is also considered a very important GPR application. However, the anticipated water table response is not observed in the GPR data at many sites (Annan et al., 1991). For example, in favourable, highly resistive weathered granite survey conditions, while working at 200 MHz GPR frequency, Mahmoudzadeh et al. (2012) could detect water table depth only down to ~
3 m b.g.s.

### 4.1.4. Seismic reflection

The seismic reflection method measures the propagation velocity of seismic waves artificially generated in rocks, soils or structures. The elastic seismic waves employ energy sources like electronically controlled vibrators, sledge hammers or controlled explosions. The velocity of propagation of initiated waves (variable for different materials crossed) is detected utilizing sensors (geophones) placed according to linear or 3D arrays on the ground or in equipped boreholes (Rubin and Hubbard, 2005). The velocity analyzed provides 3D geological and hydrostratigraphic information, including soil density and texture down to kilometres of depth, although losing resolution with depth (Bertoni et al., 2020). Active seismic using an accelerating weight striking a steel plate and geophones was shown to be able to define the position of the groundwater table (Flinchum et al., 2020). Also, employing the passive seismic, using geophones or fibre-optic cables (e.g. designed for telecommunication) showed that the measured ambient seismic noise correlates with the water table fluctuation (Garambois et al., 2019; Tribaldos and Ajo-Franklin, 2021). Like the electrical methods, the seismic reflection is affected by non-uniqueness due to the equivalence of different shear-wave velocity profiles (Foti et al., 2009).

### 4.1.5. Terrestrial gravimetry

Terrestrial gravimetry measures the acceleration of gravity at the Earth's surface, varying in space and time (Güntner et al., 2017). The major advantage of the gravity signal for hydrological studies is the direct relation between gravity variations and variations of water masses (Fores et al., 2016). Time-lapse relative gravimetry can give an integrated measure of water storage changes over tens to hundreds of cubic meters (Christiansen et al., 2011b). Extensive testing of gravimeter observation against hydrological monitoring, neutron probe measurements of water content and magnetic resonance sounding, were carried out applying absolute gravity (Hector et al., 2013) and relative gravity measurements (Pfeffer et al., 2013). The benefits of gravimeter observations for modelling total water storage changes have been assessed by Creutzfeldt et al. (2010). They calibrated a simple hydrological conceptual model against gravimeter records, soil moisture and groundwater time series and validated the model against lysimeter storage changes. They concluded that gravimeter data substantially improved predictive capability of the model. However, they also stated that the gravimeter alone could contribute neither to internal model structure nor parameterisation. Christiansen et al., (2011a) showed for the Matsibe River bank storage site at the Okavango Delta in Botswana, that time-lapse gravity data combined with hydraulic head data in a coupled hydrogeophysical inversion, allowed to decrease parameter correlation in a groundwater model and constrained specific yield, hydraulic conductivity, riverbed conductance and evapotranspiration. Also, Piccolroaz et al. (2015) showed that coupling gravimetric data with hydrological data in inverse modelling of the Alpine Vermigliana catchment, led to better constraining of the model, i.e., better identification of model parameters. The main disadvantages of the terrestrial gravimeters are high cost, significant installation effort and related lack of mobility.

### 4.1.6. Magnetic resonance sounding

MRS, sometimes also referred to as Surface Nuclear Magnetic Resonance, is the only hydrogeophysical survey that can directly provide quantitative information on the distribution of water content and permeability of rocks with depth (Boucher et al., 2009; Lubczynski and Roy, 2003; Lubczynski and Roy, 2007). However, the operational efficiency of MRS is rather low, and also the success of a survey is dependent on the geographically variable S/N ratio (Lubczynski and Roy, 2004), although that dependence has recently been mitigated by various developments in MRS hardware, software and in the ways of MRS field data acquisition (Grombacher et al., 2022; Grombacher et al., 2021;

Legchenko et al., 2016). An important improvement in field surveys is the simultaneous use of MRS with other hydrogeophysical methods, such as with ER (Skibbe et al., 2020) or with EM (Behroozmand et al., 2012; Legchenko et al., 2008; Vilhelmsen et al., 2014). Also, recent advances in signal processing by spectral analysis allowed the inversion of the subsurface with reduced model uncertainty, which permitted more reliable identification of 3D heterogeneity (Larsen et al., 2021).

In favourable S/N conditions, the current MRS hardware implementation can reach penetration depths of a maximum of 100–150 m. In low S/N conditions, MRS surveys are usually performed with discrete, 1D stratification, the penetration depth is low, and approximately only two survey locations per day can be achieved. With current MRS instrumental characteristics and relatively low time-wise operational efficiency, the MRS is not suitable for groundwater exploration. Instead, it can provide quantitative system parameterisation input for IHMs, comparable with pumping tests (Lubczynski and Roy, 2005; Vouillamoz et al., 2014) and also can be used to confirm site suitability for groundwater exploitation (Vouillamoz et al., 2002; Vouillamoz et al., 2016). At a given site, the MRS can assess porosity, specific yield, hydraulic conductivity and transmissivity, all as depth-dependent data. Recently, MRS was also used to monitor unsaturated water content (Legchenko et al., 2022) and to calibrate a HYDRUS 1D unsaturated zone model (Legchenko et al., 2020). The main disadvantage of the MRS is the high cost of the instrument.

## 4.2. Drone-borne hydrogeophysics

The last two decades have seen a rapid development of drone platforms. The technical details of drone platforms and their capabilities can be found in recent reviews by Acharya et al. (2021) and Velez-Nicolas et al. (2021). Both reviews use different terms for a flying system, the former addressing it as Unmanned Aerial Vehicle (UAV) and the latter as Unmanned Aerial System (UAS). That terminology confusion is also embedded in the scientific literature. In this study, we assume that a UAV is a drone vehicle with a navigation system, while the whole flying system, i.e., a drone with a navigation system and imaging sensors, is a UAS. The significant advantage is that it can fly on-demand equipped with modern operational sensors (e.g., GPS receiver, a WiFi serial transceiver, a distance sensor to measure the height) and various types of imaging sensors depending on the survey type. However, there are limitations for a given drone type, such as payload, power supply, wind resistance, etc. Fortunately, with rapid technological improvements, all sensors are becoming lighter while the payload, flight time, and control of drones are improving.

The six standard sensor types, typically placed on board satellites (section 3), are also being implemented onboard drones. These six are: RGB (red–green–blue), multispectral, hyperspectral, LiDAR, thermal infrared and microwave. The explanation of the advantages and disadvantages of these sensors in various drone-borne hydrological applications is discussed by Acharya et al. (2021), Velez-Nicolas et al. (2021) and Mangel et al. (2022). The most typical drone application is RGB aerial photogrammetry of the land surface, which can help identify with high precision, ground surface objects (rivers, lakes, crops, etc.) typically embedded in IHMs and monitor land cover changes. Multispectral cameras allow for determining vegetation indices, e.g. Leaf Area Index (Gong et al., 2021) or NDVI (Wahab et al., 2018), used in the assessment of rainfall interception loss (section 3.1.2) and PET (section 3.1.3), respectively. Hyperspectral sensors offer two important applications, i.e. water quality monitoring (Lu et al., 2021) and classification of plant species (Olariu et al., 2022). The main application of drone-borne thermal infrared cameras in IHMs is that they provide high-resolution surface temperature for the estimation of evapotranspiration (Niu et al., 2020), but also to identify places of surface–groundwater exchange (Abolt et al., 2018). The near-infrared LiDAR shows potential in 3D high-resolution reconstruction of surface topography and vegetation characteristics (height, roughness, and density), applicable for flood risk
3 m b.g.s.

### 4.1.4. Seismic reflection

The seismic reflection method measures the propagation velocity of seismic waves artificially generated in rocks, soils or structures. The elastic seismic waves employ energy sources like electronically controlled vibrators, sledge hammers or controlled explosions. The velocity of propagation of initiated waves (variable for different materials crossed) is detected utilizing sensors (geophones) placed according to linear or 3D arrays on the ground or in equipped boreholes (Rubin and Hubbard, 2005). The velocity analyzed provides 3D geological and hydrostratigraphic information, including soil density and texture down to kilometres of depth, although losing resolution with depth (Bertoni et al., 2020). Active seismic using an accelerating weight striking a steel plate and geophones was shown to be able to define the position of the groundwater table (Flinchum et al., 2020). Also, employing the passive seismic, using geophones or fibre-optic cables (e.g. designed for telecommunication) showed that the measured ambient seismic noise correlates with the water table fluctuation (Garambois et al., 2019; Tribaldos and Ajo-Franklin, 2021). Like the electrical methods, the seismic reflection is affected by non-uniqueness due to the equivalence of different shear-wave velocity profiles (Foti et al., 2009).

### 4.1.5. Terrestrial gravimetry

Terrestrial gravimetry measures the acceleration of gravity at the Earth's surface, varying in space and time (Güntner et al., 2017). The major advantage of the gravity signal for hydrological studies is the direct relation between gravity variations and variations of water masses (Fores et al., 2016). Time-lapse relative gravimetry can give an integrated measure of water storage changes over tens to hundreds of cubic meters (Christiansen et al., 2011b). Extensive testing of gravimeter observation against hydrological monitoring, neutron probe measurements of water content and magnetic resonance sounding, were carried out applying absolute gravity (Hector et al., 2013) and relative gravity measurements (Pfeffer et al., 2013). The benefits of gravimeter observations for modelling total water storage changes have been assessed by Creutzfeldt et al. (2010). They calibrated a simple hydrological conceptual model against gravimeter records, soil moisture and groundwater time series and validated the model against lysimeter storage changes. They concluded that gravimeter data substantially improved predictive capability of the model. However, they also stated that the gravimeter alone could contribute neither to internal model structure nor parameterisation. Christiansen et al., (2011a) showed for the Matsibe River bank storage site at the Okavango Delta in Botswana, that time-lapse gravity data combined with hydraulic head data in a coupled hydrogeophysical inversion, allowed to decrease parameter correlation in a groundwater model and constrained specific yield, hydraulic conductivity, riverbed conductance and evapotranspiration. Also, Piccolroaz et al. (2015) showed that coupling gravimetric data with hydrological data in inverse modelling of the Alpine Vermigliana catchment, led to better constraining of the model, i.e., better identification of model parameters. The main disadvantages of the terrestrial gravimeters are high cost, significant installation effort and related lack of mobility.

### 4.1.6. Magnetic resonance sounding

MRS, sometimes also referred to as Surface Nuclear Magnetic Resonance, is the only hydrogeophysical survey that can directly provide quantitative information on the distribution of water content and permeability of rocks with depth (Boucher et al., 2009; Lubczynski and Roy, 2003; Lubczynski and Roy, 2007). However, the operational efficiency of MRS is rather low, and also the success of a survey is dependent on the geographically variable S/N ratio (Lubczynski and Roy, 2004), although that dependence has recently been mitigated by various developments in MRS hardware, software and in the ways of MRS field data acquisition (Grombacher et al., 2022; Grombacher et al., 2021;

Legchenko et al., 2016). An important improvement in field surveys is the simultaneous use of MRS with other hydrogeophysical methods, such as with ER (Skibbe et al., 2020) or with EM (Behroozmand et al., 2012; Legchenko et al., 2008; Vilhelmsen et al., 2014). Also, recent advances in signal processing by spectral analysis allowed the inversion of the subsurface with reduced model uncertainty, which permitted more reliable identification of 3D heterogeneity (Larsen et al., 2021).

In favourable S/N conditions, the current MRS hardware implementation can reach penetration depths of a maximum of 100–150 m. In low S/N conditions, MRS surveys are usually performed with discrete, 1D stratification, the penetration depth is low, and approximately only two survey locations per day can be achieved. With current MRS instrumental characteristics and relatively low time-wise operational efficiency, the MRS is not suitable for groundwater exploration. Instead, it can provide quantitative system parameterisation input for IHMs, comparable with pumping tests (Lubczynski and Roy, 2005; Vouillamoz et al., 2014) and also can be used to confirm site suitability for groundwater exploitation (Vouillamoz et al., 2002; Vouillamoz et al., 2016). At a given site, the MRS can assess porosity, specific yield, hydraulic conductivity and transmissivity, all as depth-dependent data. Recently, MRS was also used to monitor unsaturated water content (Legchenko et al., 2022) and to calibrate a HYDRUS 1D unsaturated zone model (Legchenko et al., 2020). The main disadvantage of the MRS is the high cost of the instrument.

## 4.2. Drone-borne hydrogeophysics

The last two decades have seen a rapid development of drone platforms. The technical details of drone platforms and their capabilities can be found in recent reviews by Acharya et al. (2021) and Velez-Nicolas et al. (2021). Both reviews use different terms for a flying system, the former addressing it as Unmanned Aerial Vehicle (UAV) and the latter as Unmanned Aerial System (UAS). That terminology confusion is also embedded in the scientific literature. In this study, we assume that a UAV is a drone vehicle with a navigation system, while the whole flying system, i.e., a drone with a navigation system and imaging sensors, is a UAS. The significant advantage is that it can fly on-demand equipped with modern operational sensors (e.g., GPS receiver, a WiFi serial transceiver, a distance sensor to measure the height) and various types of imaging sensors depending on the survey type. However, there are limitations for a given drone type, such as payload, power supply, wind resistance, etc. Fortunately, with rapid technological improvements, all sensors are becoming lighter while the payload, flight time, and control of drones are improving.

The six standard sensor types, typically placed on board satellites (section 3), are also being implemented onboard drones. These six are: RGB (red–green–blue), multispectral, hyperspectral, LiDAR, thermal infrared and microwave. The explanation of the advantages and disadvantages of these sensors in various drone-borne hydrological applications is discussed by Acharya et al. (2021), Velez-Nicolas et al. (2021) and Mangel et al. (2022). The most typical drone application is RGB aerial photogrammetry of the land surface, which can help identify with high precision, ground surface objects (rivers, lakes, crops, etc.) typically embedded in IHMs and monitor land cover changes. Multispectral cameras allow for determining vegetation indices, e.g. Leaf Area Index (Gong et al., 2021) or NDVI (Wahab et al., 2018), used in the assessment of rainfall interception loss (section 3.1.2) and PET (section 3.1.3), respectively. Hyperspectral sensors offer two important applications, i.e. water quality monitoring (Lu et al., 2021) and classification of plant species (Olariu et al., 2022). The main application of drone-borne thermal infrared cameras in IHMs is that they provide high-resolution surface temperature for the estimation of evapotranspiration (Niu et al., 2020), but also to identify places of surface–groundwater exchange (Abolt et al., 2018). The near-infrared LiDAR shows potential in 3D high-resolution reconstruction of surface topography and vegetation characteristics (height, roughness, and density), applicable for flood risk
assessment and for monitoring ecosystem changes (Resop et al., 2021; Trepekli et al., 2022). However, it has been shown (Bandini et al., 2020) that its use for the assessment of water bodies is not optimal, being outperformed by green (bathymetric) LiDAR (Mano et al., 2020; Szafarczyk and Toś, 2023) and GPR. Finally, there are implementations of passive L-band microwave sensors (in contrast to active GPR), mounted on board of drones to monitor soil moisture (Acevo-Herrera et al., 2010; Houtz et al., 2020a; Houtz et al., 2020b).

Recent technological advances are observed in adopting classic surface hydrogeophysical methods onboard drones. The most significant such drone developments are with GPR and its applications in soil moisture mapping. In that respect, Wu et al. (2019) and Wu et al. (2022) proposed a whole radar system with a lightweight antenna of only 1.5 kg, based on a handheld vector network analyzer, working in the frequency domain and applying a full-wave inversion to derive soil moisture. Another drone GPR application is presented by Valence et al. (2022), who conducted a snowpack evolution survey in a 200 m² area in Canada, with weekly surveys of snow depth, density, snow water equivalent and bulk liquid water content. Compared to surface GPR, the advantages of such a drone-borne survey are that the snow is not compacted, and the survey is faster. Finally, a very interesting, drone-borne GPR application for mapping inland water bathymetry is presented by Bandini et al. (2023). Their study confirmed larger efficiency and similar accuracy as a GPR survey from a boat (see 4.1.3), but had two limitations: i) a more restrictive minimum water body depth requirement (typically 0.8–1.1 m for drone-borne GPR, while 0.3–0.4 m for boat GPR); ii) requirement to fly the GPR antenna at altitudes of approx. 0.5 m above the water surface. Not only GPR but also the electromagnetic method has its recent drone-borne implementations. The version developed in The Netherlands (Karaoulis et al., 2022) focused on low-cost monitoring of local water saturation and salinity changes, while the version developed in Japan (Mitsuhata et al., 2022) aimed at soil resistivity mapping.

## 4.3. Airborne hydrogeophysics

There are hydrogeophysical surveys that cannot be performed with drones but only with piloted helicopters or airplanes, especially when: i) heavy instruments need to be onboard (large payload required); ii) a large area needs to be flown (aircraft can fly at high altitude and have sufficient power supply); iii) harsh conditions of survey prevail (aircraft can overcome significant wind or fly in inaccessible areas); iv) survey need to be done fast (better speed and maneuverability).

The same six standard sensor types as used onboard drones and discussed in the previous section 4.2, have also been widely utilised in airborne campaigns. Airborne surveys with these sensors are similar to drone-borne surveys, differing mainly by larger airborne-surveyed areas; therefore, their descriptions are not repeated here. However, in contrast to drones introduced only recently, airborne campaigns with these six standard sensors (or at least some of them) have been carried out since the early 20th Century, so some archive data is available in many regions. For example, in Niger, historical panchromatic aerial photographs enabled change detection from 1950 to 1992 and showed the important impact of deforestation on both, surface runoff and groundwater recharge (Leblanc et al., 2008).

Considering classic hydrogeophysical airborne applications, the airborne EM method is the most frequently used for groundwater applications and, hence, most applicable for IHMs. The airborne EM can operate as FDEM or TDEM (see 4.1.2) from helicopters or airplanes equipped with a transmitter loop or coil to induce current flow into the ground and create a secondary magnetic field that a receiver coil senses. The airborne EM signals can reach depths of several hundred meters, depending on the method used, flight height, type of instrument, and ground conductivity (Christiansen et al., 2006; Siemon, 2009). The main applications of airborne EM include: i) water quality mapping, in particular, to delineate fresh and saltwater zones in aquifers (Chongo et al., 2015; Podgorski et al., 2017); ii) acquisition of lithological information applied for 3D subsurface hydrostratigraphy (Marker et al., 2015); iii) determination of 3D spatial hydraulic conductivity distribution, e.g. at the catchment scale, integrating ground measurements of electrical resistivity and borehole data (O'Connell et al., 2020); iv) identification of groundwater resources, e.g. in buried valleys (Morgan et al., 2019).

The airborne EM campaigns are typically combined with simultaneous magnetic and gamma-ray measurements, which in contrast to airborne EM, do not directly contribute to groundwater assessments, but still can provide useful data for IHMs. The magnetometer measurement can delineate igneous rock intrusions at depth (e.g. lava flows such as basalt), local magnetic anomalies caused by faults or dykes, etc. Gamma-ray spectrometers measure radioelement concentrations in a thin superficial layer of the earth's surface. It is applied mainly for geology, lithology, regolith and soil mapping (Rubin and Hubbard, 2005). The SkyTEM system is a popular example of an airborne implementation of TEM, magnetometer and gamma-ray spectrometer, all mounted directly on the TEM frame. The SkyTEM system is typically carried by a helicopter, allowing measurements close to the ground surface (Steuer et al., 2009).

There are also airborne GPR implementations. Easily accessible areas can be operated by drone GPR systems, but difficult-to-access areas, like mountainous areas, would need to be operated by helicopter-borne GPR. Most of such applications are related to the estimation of ice thickness and below-ice morphology. Pros and cons of helicopter-borne GPR in rugged mountainous areas are discussed by (Forte et al., 2019).

## 4.4. Integrations of hydrogeophysics with models

Various hydrogeophysical techniques in support of conceptual and numerical groundwater models have been well documented (Dam and Christensen, 2003). Considering surface hydrogeophysics, for example, GPR was used to develop a conceptual model of an aquifer architecture (facies heterogeneity) within a shallow coastal alluvial plain (Ezzy et al., 2006). In contrast, TDEM was used to calibrate a saltwater intrusion model (Herckenrath et al., 2013b) and reduce groundwater model prediction error (Christensen et al., 2016). Also, airborne hydrogeophysics contributed to groundwater models. Dickson et al. (2014) utilised airborne and ground-based magnetic surveys to map low-permeability dolerite dykes for implementation in a groundwater flow model. Morgan et al. (2019) used airborne electromagnetic (SkyTEM), borehole gamma-ray and lithology logs to map buried valleys and conceptualize a groundwater flow model.

Many groundwater modelling studies used two or more hydrogeophysical methods simultaneously. For example, Herckenrath et al., (2013a) applied sequential and joint hydrogeophysical inversion of TDEM and ER tomography in a groundwater flow model. Boucher et al. (2012) used MRS supported by TDEM to constrain model calibration and reduce model uncertainty. Luoma et al. (2021) conducted a complex study applying many different methods such as gravimetric, ground-penetrating radar, shallow seismic surveys, drill logs, groundwater level monitoring data, and a LiDAR digital elevation model to setup a groundwater flow model to quantify submarine groundwater discharge.

The use of hydrogeophysics for IHMs is relatively new. Baroncini-Turricchia et al. (2014) used MRS to support the hydrogeological parameterisation of a MARMITES-MODFLOW IHM of the Carrizal catchment, Spain. Marker et al. (2015) defined the hydrostratigraphy of a MIKE SHE IHM using airborne EM and lithological borehole logs. Ala-Aho et al. (2015) used airborne thermal infrared imaging to define the spatial occurrence of groundwater-lake interactions, implemented in the HydroGeoSphere IHM. Lesparre et al. (2020) conditioned the Normally Integrated Hydrological Model with MRS measurements and streamflow data to provide reliable model calibration in a data-scarce catchment. Various geophysical techniques, including MRS, GPR, FDEM and ER, next to RS and field observations, were applied by
assessment and for monitoring ecosystem changes (Resop et al., 2021; Trepekli et al., 2022). However, it has been shown (Bandini et al., 2020) that its use for the assessment of water bodies is not optimal, being outperformed by green (bathymetric) LiDAR (Mano et al., 2020; Szafarczyk and Toś, 2023) and GPR. Finally, there are implementations of passive L-band microwave sensors (in contrast to active GPR), mounted on board of drones to monitor soil moisture (Acevo-Herrera et al., 2010; Houtz et al., 2020a; Houtz et al., 2020b).

Recent technological advances are observed in adopting classic surface hydrogeophysical methods onboard drones. The most significant such drone developments are with GPR and its applications in soil moisture mapping. In that respect, Wu et al. (2019) and Wu et al. (2022) proposed a whole radar system with a lightweight antenna of only 1.5 kg, based on a handheld vector network analyzer, working in the frequency domain and applying a full-wave inversion to derive soil moisture. Another drone GPR application is presented by Valence et al. (2022), who conducted a snowpack evolution survey in a 200 m² area in Canada, with weekly surveys of snow depth, density, snow water equivalent and bulk liquid water content. Compared to surface GPR, the advantages of such a drone-borne survey are that the snow is not compacted, and the survey is faster. Finally, a very interesting, drone-borne GPR application for mapping inland water bathymetry is presented by Bandini et al. (2023). Their study confirmed larger efficiency and similar accuracy as a GPR survey from a boat (see 4.1.3), but had two limitations: i) a more restrictive minimum water body depth requirement (typically 0.8–1.1 m for drone-borne GPR, while 0.3–0.4 m for boat GPR); ii) requirement to fly the GPR antenna at altitudes of approx. 0.5 m above the water surface. Not only GPR but also the electromagnetic method has its recent drone-borne implementations. The version developed in The Netherlands (Karaoulis et al., 2022) focused on low-cost monitoring of local water saturation and salinity changes, while the version developed in Japan (Mitsuhata et al., 2022) aimed at soil resistivity mapping.

## 4.3. Airborne hydrogeophysics

There are hydrogeophysical surveys that cannot be performed with drones but only with piloted helicopters or airplanes, especially when: i) heavy instruments need to be onboard (large payload required); ii) a large area needs to be flown (aircraft can fly at high altitude and have sufficient power supply); iii) harsh conditions of survey prevail (aircraft can overcome significant wind or fly in inaccessible areas); iv) survey need to be done fast (better speed and maneuverability).

The same six standard sensor types as used onboard drones and discussed in the previous section 4.2, have also been widely utilised in airborne campaigns. Airborne surveys with these sensors are similar to drone-borne surveys, differing mainly by larger airborne-surveyed areas; therefore, their descriptions are not repeated here. However, in contrast to drones introduced only recently, airborne campaigns with these six standard sensors (or at least some of them) have been carried out since the early 20th Century, so some archive data is available in many regions. For example, in Niger, historical panchromatic aerial photographs enabled change detection from 1950 to 1992 and showed the important impact of deforestation on both, surface runoff and groundwater recharge (Leblanc et al., 2008).

Considering classic hydrogeophysical airborne applications, the airborne EM method is the most frequently used for groundwater applications and, hence, most applicable for IHMs. The airborne EM can operate as FDEM or TDEM (see 4.1.2) from helicopters or airplanes equipped with a transmitter loop or coil to induce current flow into the ground and create a secondary magnetic field that a receiver coil senses. The airborne EM signals can reach depths of several hundred meters, depending on the method used, flight height, type of instrument, and ground conductivity (Christiansen et al., 2006; Siemon, 2009). The main applications of airborne EM include: i) water quality mapping, in particular, to delineate fresh and saltwater zones in aquifers (Chongo et al., 2015; Podgorski et al., 2017); ii) acquisition of lithological information applied for 3D subsurface hydrostratigraphy (Marker et al., 2015); iii) determination of 3D spatial hydraulic conductivity distribution, e.g. at the catchment scale, integrating ground measurements of electrical resistivity and borehole data (O'Connell et al., 2020); iv) identification of groundwater resources, e.g. in buried valleys (Morgan et al., 2019).

The airborne EM campaigns are typically combined with simultaneous magnetic and gamma-ray measurements, which in contrast to airborne EM, do not directly contribute to groundwater assessments, but still can provide useful data for IHMs. The magnetometer measurement can delineate igneous rock intrusions at depth (e.g. lava flows such as basalt), local magnetic anomalies caused by faults or dykes, etc. Gamma-ray spectrometers measure radioelement concentrations in a thin superficial layer of the earth's surface. It is applied mainly for geology, lithology, regolith and soil mapping (Rubin and Hubbard, 2005). The SkyTEM system is a popular example of an airborne implementation of TEM, magnetometer and gamma-ray spectrometer, all mounted directly on the TEM frame. The SkyTEM system is typically carried by a helicopter, allowing measurements close to the ground surface (Steuer et al., 2009).

There are also airborne GPR implementations. Easily accessible areas can be operated by drone GPR systems, but difficult-to-access areas, like mountainous areas, would need to be operated by helicopter-borne GPR. Most of such applications are related to the estimation of ice thickness and below-ice morphology. Pros and cons of helicopter-borne GPR in rugged mountainous areas are discussed by (Forte et al., 2019).

## 4.4. Integrations of hydrogeophysics with models

Various hydrogeophysical techniques in support of conceptual and numerical groundwater models have been well documented (Dam and Christensen, 2003). Considering surface hydrogeophysics, for example, GPR was used to develop a conceptual model of an aquifer architecture (facies heterogeneity) within a shallow coastal alluvial plain (Ezzy et al., 2006). In contrast, TDEM was used to calibrate a saltwater intrusion model (Herckenrath et al., 2013b) and reduce groundwater model prediction error (Christensen et al., 2016). Also, airborne hydrogeophysics contributed to groundwater models. Dickson et al. (2014) utilised airborne and ground-based magnetic surveys to map low-permeability dolerite dykes for implementation in a groundwater flow model. Morgan et al. (2019) used airborne electromagnetic (SkyTEM), borehole gamma-ray and lithology logs to map buried valleys and conceptualize a groundwater flow model.

Many groundwater modelling studies used two or more hydrogeophysical methods simultaneously. For example, Herckenrath et al., (2013a) applied sequential and joint hydrogeophysical inversion of TDEM and ER tomography in a groundwater flow model. Boucher et al. (2012) used MRS supported by TDEM to constrain model calibration and reduce model uncertainty. Luoma et al. (2021) conducted a complex study applying many different methods such as gravimetric, ground-penetrating radar, shallow seismic surveys, drill logs, groundwater level monitoring data, and a LiDAR digital elevation model to setup a groundwater flow model to quantify submarine groundwater discharge.

The use of hydrogeophysics for IHMs is relatively new. Baroncini-Turricchia et al. (2014) used MRS to support the hydrogeological parameterisation of a MARMITES-MODFLOW IHM of the Carrizal catchment, Spain. Marker et al. (2015) defined the hydrostratigraphy of a MIKE SHE IHM using airborne EM and lithological borehole logs. Ala-Aho et al. (2015) used airborne thermal infrared imaging to define the spatial occurrence of groundwater-lake interactions, implemented in the HydroGeoSphere IHM. Lesparre et al. (2020) conditioned the Normally Integrated Hydrological Model with MRS measurements and streamflow data to provide reliable model calibration in a data-scarce catchment. Various geophysical techniques, including MRS, GPR, FDEM and ER, next to RS and field observations, were applied by
Francés et al. (2014) to define a conceptual model of an experimental, about 80 km² granitic Sardon catchment, recently implemented in the MODFLOW 6 IHM by Daoud et al. (2022).

## 5. Discussion

Traditionally, the predictive capacities of standalone (single hydrological domain) groundwater, unsaturated zone or land-surface models have always suffered from large uncertainty, mainly due to limited knowledge of boundary conditions. Standalone models used to have particularly poor performance in shallow groundwater table environments with distinct surface-groundwater interactions, where frequent and important water flux exchanges between aquifer, unsaturated zone, land surface, and the atmosphere occur. Still, two decades back, Bateelaan and De Smedt (2004) proposed groundwater model improvements in that respect, but even these improvements could not quantify surface-groundwater interactions adequately because standalone models are not adapted to such tasks.

The complexity of water flux exchanges across the unsaturated zone, especially at its bottom (net recharge) and top (infiltration and exfiltration), stimulated the replacement of standalone models by IHMs. That integration of surface water domain with groundwater domain across the unsaturated zone makes the IHM setup and calibration more complex but also more reliable because: i) water fluxes of all model domains are dynamically integrated into one hydrological system; ii) the top IHM boundary, with its spatiotemporally variable climate forcings, can be defined in a more realistic manner than the top net recharge boundary in standalone groundwater models, especially in shallow water table environments; iii) the top and the bottom of the IHM unsaturated zone domain are represented more realistically, e.g. by field, RS and hydrogeophysical measurements, than the bottom in standalone land-surface and the top in groundwater models; iv) many different types of observations, also spatiotemporal RS and hydrogeophysical data covering the entire grid domain, can be used to constrain IHM state variables, in a multivariate calibration or data assimilation.

RS is a rapidly expanding field of science providing space-borne data at continuously improving spatiotemporal resolution. Satellite products (e.g., precipitation, PET, ET, SM, stages etc.) have become widely available online. Their use typically does not require specialised RS-processing knowledge, so RS products can be directly downloaded from online data portals. The most widely used platforms to access archive remote sensing data include NASA (https://www.earthdata.nasa.gov/), ESA (https://earth.esa.int/eogateway/) and Google Earth Engine (https://earthengine.google.com). Moreover, a user can also have access to near-real time RS products, for example, Copernicus sites (https://land.copernicus.eu/ and https://dataspace.copernicus.eu/), but also through self-installed, cost-effective satellite-based dissemination systems such as GEONETCast, which facilitates easy import of access-free and open ftp/http time series of environmental data (Maathuis et al., 2012).

The main advantages of RS products as input for IHMs are the spatially and temporally semi-continuous data, with historical archives (change detection), typically free-of-charge. The main disadvantages are uncertainty and frequently insufficient spatial resolution compared to the heterogeneity and complexity of an area simulated by an IHM. However, the uncertainty of RS products can be substantially reduced by their integration (mainly bias-correction) with in-situ data (Gebremedhin et al., 2021) or, for some applications, also with drone-borne measurements.

The main advantages of airborne and drone-borne hydrogeophysics as compared to RS are that they can: i) observe not only the ground surface but also the subsurface, even deep subsurface; ii) provide high spatial and temporal resolution typically higher than RS; iii) acquire data 'on demand' according to project need, including in cloudy conditions; and iv) easily integrate with ground measurements. The main disadvantages of airborne and drone-borne hydrogeophysics as

compared to RS are: i) the higher costs of campaigns, although less for those using drones; ii) lower spatial efficiency with a smaller footprint; iii) requirement of specialised knowledge; iv) requirement of cumbersome data processing; v) lack of long time series archives. Compared to surface hydrogeophysics, the main advantages of airborne and drone-borne hydrogeophysics are their better efficiency and larger footprint. The main disadvantages are the lower spatial resolution and, hence, the lower ability to resolve small targets. Also, not all the hydrogeophysical methods are adapted to airborne and drone-borne surveys, e.g. MRS is not.

IHMs use four types of input data (Fig. 1): i) climate forcings (spatiotemporally variable); ii) parameters (spatially variable); iii) boundary conditions (spatially or spatiotemporally variable); and iv) observations to constrain state variables (spatiotemporally variable). Climate forcings, parameters, and boundary conditions are indispensable to running an IHM, which produces simulated state variables that are compared with observations to minimize the differences between them. Climate forcings and some boundary conditions are typically required daily as continuous time series, in contrast to state variables that can be assigned incidentally. Climate forcings must be assigned over the entire model domain, parameters in every grid cell at any layer, while boundary conditions and state variables can be set locally. These requirements constrain different data sources and methods of IHM input data preprocessing.

### 5.1. Climate forcings

Precipitation is the most important input data of IHMs and the most difficult to define properly because of its large spatiotemporal variability. As transient IHMs require time series of precipitation covering the entire grid and every stress period simulated, a combination of satellite RS products and ground measurements is often used for that purpose. Inaccurate precipitation estimates, e.g. made by RS only, can jeopardize IHMs as the precipitation uncertainty propagates into the IHM results. Besides, standard 0.1° RS precipitation products (or coarser) have inherent uncertainty related to their large pixel area (typically ≥ 100 km²), particularly distinct in topographically complex terrains. In such pixels, the comparison of the mean precipitation with the point (gauge) precipitation estimate is problematic (Jiang and Bauer-Gottwein, 2019). The true precipitation at coarse scale products (10–100 km) is essentially unknown. Therefore, higher spatiotemporal resolution products (e.g. MPEG with ~ 3 km spatial and 15 min temporal resolution), especially after bias correction using ground measurements, in general, provide more reliable results.

Numerous RS precipitation products perform differently in various locations, as presented in section 3.1.1. Therefore prior to using a particular RS precipitation product in an IHM, it is recommended to: i) validate different RS precipitation products against available rain gauges to select the most reliable RS product; if there are no rain gauges or an insufficient number to provide reliable evaluation, one may consider installing some, even for a limited period, to obtain information, just for validation and bias correction; ii) select and apply an optimal bias correction method accounting for the complexity of an investigated area; iii) validate the bias-corrected RS precipitation; iv) if necessary scale (down or up) the RS precipitation to match it with the IHM grid.

Rainfall interception loss can represent a large part of precipitation but is typically underestimated or even neglected by modellers. If interception represents even only a few percent of rainfall, it can substantially change the water budget of a model (Hassan and Lubczynski, 2024). However, in forests of tropical and moderate climates, interception can exceed even 30 % of rainfall (Gerrits et al., 2010; McJannet et al., 2007; Miralles et al., 2010), while in seasonally dry savannah, it largely depends on tree density and species composition (evergreen or deciduous), but together with the forest floor it can also exceed 30 % (Tsiko et al., 2012). So, interception is an essential component of IHMs
Francés et al. (2014) to define a conceptual model of an experimental, about 80 km² granitic Sardon catchment, recently implemented in the MODFLOW 6 IHM by Daoud et al. (2022).

## 5. Discussion

Traditionally, the predictive capacities of standalone (single hydrological domain) groundwater, unsaturated zone or land-surface models have always suffered from large uncertainty, mainly due to limited knowledge of boundary conditions. Standalone models used to have particularly poor performance in shallow groundwater table environments with distinct surface-groundwater interactions, where frequent and important water flux exchanges between aquifer, unsaturated zone, land surface, and the atmosphere occur. Still, two decades back, Bateelaan and De Smedt (2004) proposed groundwater model improvements in that respect, but even these improvements could not quantify surface-groundwater interactions adequately because standalone models are not adapted to such tasks.

The complexity of water flux exchanges across the unsaturated zone, especially at its bottom (net recharge) and top (infiltration and exfiltration), stimulated the replacement of standalone models by IHMs. That integration of surface water domain with groundwater domain across the unsaturated zone makes the IHM setup and calibration more complex but also more reliable because: i) water fluxes of all model domains are dynamically integrated into one hydrological system; ii) the top IHM boundary, with its spatiotemporally variable climate forcings, can be defined in a more realistic manner than the top net recharge boundary in standalone groundwater models, especially in shallow water table environments; iii) the top and the bottom of the IHM unsaturated zone domain are represented more realistically, e.g. by field, RS and hydrogeophysical measurements, than the bottom in standalone land-surface and the top in groundwater models; iv) many different types of observations, also spatiotemporal RS and hydrogeophysical data covering the entire grid domain, can be used to constrain IHM state variables, in a multivariate calibration or data assimilation.

RS is a rapidly expanding field of science providing space-borne data at continuously improving spatiotemporal resolution. Satellite products (e.g., precipitation, PET, ET, SM, stages etc.) have become widely available online. Their use typically does not require specialised RS-processing knowledge, so RS products can be directly downloaded from online data portals. The most widely used platforms to access archive remote sensing data include NASA (https://www.earthdata.nasa.gov/), ESA (https://earth.esa.int/eogateway/) and Google Earth Engine (https://earthengine.google.com). Moreover, a user can also have access to near-real time RS products, for example, Copernicus sites (https://land.copernicus.eu/ and https://dataspace.copernicus.eu/), but also through self-installed, cost-effective satellite-based dissemination systems such as GEONETCast, which facilitates easy import of access-free and open ftp/http time series of environmental data (Maathuis et al., 2012).

The main advantages of RS products as input for IHMs are the spatially and temporally semi-continuous data, with historical archives (change detection), typically free-of-charge. The main disadvantages are uncertainty and frequently insufficient spatial resolution compared to the heterogeneity and complexity of an area simulated by an IHM. However, the uncertainty of RS products can be substantially reduced by their integration (mainly bias-correction) with in-situ data (Gebremedhin et al., 2021) or, for some applications, also with drone-borne measurements.

The main advantages of airborne and drone-borne hydrogeophysics as compared to RS are that they can: i) observe not only the ground surface but also the subsurface, even deep subsurface; ii) provide high spatial and temporal resolution typically higher than RS; iii) acquire data 'on demand' according to project need, including in cloudy conditions; and iv) easily integrate with ground measurements. The main disadvantages of airborne and drone-borne hydrogeophysics as

compared to RS are: i) the higher costs of campaigns, although less for those using drones; ii) lower spatial efficiency with a smaller footprint; iii) requirement of specialised knowledge; iv) requirement of cumbersome data processing; v) lack of long time series archives. Compared to surface hydrogeophysics, the main advantages of airborne and drone-borne hydrogeophysics are their better efficiency and larger footprint. The main disadvantages are the lower spatial resolution and, hence, the lower ability to resolve small targets. Also, not all the hydrogeophysical methods are adapted to airborne and drone-borne surveys, e.g. MRS is not.

IHMs use four types of input data (Fig. 1): i) climate forcings (spatiotemporally variable); ii) parameters (spatially variable); iii) boundary conditions (spatially or spatiotemporally variable); and iv) observations to constrain state variables (spatiotemporally variable). Climate forcings, parameters, and boundary conditions are indispensable to running an IHM, which produces simulated state variables that are compared with observations to minimize the differences between them. Climate forcings and some boundary conditions are typically required daily as continuous time series, in contrast to state variables that can be assigned incidentally. Climate forcings must be assigned over the entire model domain, parameters in every grid cell at any layer, while boundary conditions and state variables can be set locally. These requirements constrain different data sources and methods of IHM input data preprocessing.

### 5.1. Climate forcings

Precipitation is the most important input data of IHMs and the most difficult to define properly because of its large spatiotemporal variability. As transient IHMs require time series of precipitation covering the entire grid and every stress period simulated, a combination of satellite RS products and ground measurements is often used for that purpose. Inaccurate precipitation estimates, e.g. made by RS only, can jeopardize IHMs as the precipitation uncertainty propagates into the IHM results. Besides, standard 0.1° RS precipitation products (or coarser) have inherent uncertainty related to their large pixel area (typically ≥ 100 km²), particularly distinct in topographically complex terrains. In such pixels, the comparison of the mean precipitation with the point (gauge) precipitation estimate is problematic (Jiang and Bauer-Gottwein, 2019). The true precipitation at coarse scale products (10–100 km) is essentially unknown. Therefore, higher spatiotemporal resolution products (e.g. MPEG with ~ 3 km spatial and 15 min temporal resolution), especially after bias correction using ground measurements, in general, provide more reliable results.

Numerous RS precipitation products perform differently in various locations, as presented in section 3.1.1. Therefore prior to using a particular RS precipitation product in an IHM, it is recommended to: i) validate different RS precipitation products against available rain gauges to select the most reliable RS product; if there are no rain gauges or an insufficient number to provide reliable evaluation, one may consider installing some, even for a limited period, to obtain information, just for validation and bias correction; ii) select and apply an optimal bias correction method accounting for the complexity of an investigated area; iii) validate the bias-corrected RS precipitation; iv) if necessary scale (down or up) the RS precipitation to match it with the IHM grid.

Rainfall interception loss can represent a large part of precipitation but is typically underestimated or even neglected by modellers. If interception represents even only a few percent of rainfall, it can substantially change the water budget of a model (Hassan and Lubczynski, 2024). However, in forests of tropical and moderate climates, interception can exceed even 30 % of rainfall (Gerrits et al., 2010; McJannet et al., 2007; Miralles et al., 2010), while in seasonally dry savannah, it largely depends on tree density and species composition (evergreen or deciduous), but together with the forest floor it can also exceed 30 % (Tsiko et al., 2012). So, interception is an essential component of IHMs
and must not be disregarded in IHMs, particularly as there are nowadays RS and hydrogeophysical tools to estimate it even in environments without in-situ measurements. Such estimates are by far more reliable than arbitrary interception assumptions, especially those neglecting interception.

Field measurements of interception (see 3.1.2) are tedious and require lots of gauges, especially in areas with complex vegetation patterns and composition; as such, they are accurate only at the plot scale. At the catchment (or coarser) scale, if interception measurements per each plant species are available, spatiotemporal interception can be scaled up by applying very high-resolution RS (Hassan et al., 2017), but also by airborne or drone-borne hydrogeophysics. However, if field interception measurements are not available or available only for some species, a use of RS (Cui and Jia, 2014; Cui et al., 2015; Gebremedhin et al., 2023) or drone-borne (Bolanos-Sánchez et al., 2021) driven analytical interception model, such as the Gash model, is recommended.

Most IHMs require as input spatiotemporally variable PET forcing (see 3.1.3). Although the PET definition is quite well-accepted and used, there is no consensus on the method for its estimation. For the PET, after Dassargues (2018), we recommend PET = ETo * Kc following the FAO approach because: i) the FAO ETo is standardised (Allen et al., 1998) and can be well defined not only by ground measurements but also by RS estimates (Gebremedhin et al., 2022); ii) the Kc for different crops is already well defined (Allen et al., 1998; Pereira et al., 2021a; Pereira et al., 2021b); iii) the Kc of natural vegetation can be estimated from spatiotemporally variable indices, such as e.g., from the NDVI, which nowadays can be derived not only from RS and airborne hydrogeophysics but also at very high resolution using Unmanned Aerial Systems (Jiang et al., 2020).

## 5.2. Parameters

IHMs require parameters of all hydrological domains. Satellite RS provides land surface and shallow soil characterisation but is of little or no use for subsurface hydrostratigraphy and parameterisation. Optical RS imagery provides useful information for setting up an IHM, such as catchment boundaries, slopes, rivers and water bodies' locations, etc. while multi-temporal optical data allows monitoring of important land surface changes (land use-land cover change, temporal variability of vegetation indices, etc.). Surface topography (DEM) is a very important IHM contribution of RS (section 3.2), drone-borne (section 4.2) and airborne (section 4.3) hydrogeophysics. However, one aspect of DEM is often underrated in setting and calibrating hydrological models: elevation accuracy. The accuracy of elevation greatly influences the reliability of IHMs because the estimates of hydraulic heads and stages of surface water bodies, directly depend on the precision of the elevation estimates. For example, how precise is a hydraulic head calculated from field measurements of water table depth with +/- 1 cm accuracy and RS elevation estimate with +/- ~10 m? Hence, standard RS DEM products with elevation accuracy of several meters must not be used for referencing hydraulic heads and stages of surface water bodies. The best elevation accuracy can obviously be obtained on the ground with geodetic, differential GPS providing accuracy of +/-1 cm, or by cumbersome, but also precise, classic geodmetrical equipment. A good option is also airborne or drone-borne green LiDAR or radar altimetry (see 4.2 and 4.3), especially in areas not accessible for ground-based elevation measurements. Considering satellite data, to our knowledge, the only comparable elevation accuracy is provided by the data of the ICESat-2 satellite mission, but its data availability declines from the poles towards the equator (see 3.3).

Parameterisation of shallow soils refers to soil type and texture, which often leads to the development of pedo-transfer functions. The microwave ability of soil penetration has been widely used for that purpose. Currently, the microwave Sentinel-1 operating in C-band (5.4 GHz) is the most suitable for that purpose, but in 2024, very high spatial resolution NISAR (NASA-ISRO Synthetic Aperture Radar) will be

launched, which will acquire data in S-band (3.2 GHz) and L-band (1.25 GHz). Soil type and texture can also be adequately provided by airborne and drone-borne hydrogeophysical platforms (section 4.2 and 4.3). In that respect, excellent results have been obtained from multi-temporal airborne microwave data acquired at 1.3 GHz and in different polarisations by Marzahn and Meyer (2020) as validated at the TERENO test-site "North-Eastern German Lowland Observatory", resulting in a mean RMSE of only 2.42 (Mass-%). The main disadvantage of microwave soil parameterisation is the limited signal penetration, which even for L-band is in order of 3–5 cm (Xu et al., 2014).

Parameterisation of the medium and deep subsurface is optimally done by combining hydrogeophysical and borehole data, while RS contribution is marginal. The hydraulic properties of the unsaturated zone can be investigated by surface hydrogeophysics (Francés and Lubczynski, 2011; Gallistl et al., 2022), but also more efficiently by airborne and recently also by drone-borne EM and GPR methods, which can also be used for spatial extrapolation of borehole data and surface hydrogeophysical surveys. Considering quantitative hydraulic parameterisation of the subsurface (unsaturated and saturated zone), the surface hydrogeophysical method of MRS is the only method providing quantitative estimates of hydrogeological parameters such as aquifer transmissivity, hydraulic conductivity and aquifer storage coefficient.

Hydrostratigraphy is an important element of conceptualised 3D pictures of the subsurface (Enemark et al., 2019). As such, it is indispensable in all constructed IHMs. 3D hydrostratigraphic models involving hydrostratigraphy and parameterisation, can be built based on borehole logs if they are sufficiently available. However, in many projects, the density of borehole logs is insufficient. Hydrogeophysical methods can provide supplementary hydrostratigraphic information to fill data gaps between the boreholes, although with different contributions, accuracy, and efficiency in different media (Enemark et al., 2020). The EM, ER, GPR and seismic reflection methods are affected by non-uniqueness mainly because the relationships between their measured physical signal and hydrogeological interpretation are proxies, confusing different hydrogeological media of similar physical properties, but they are operationally efficient. Therefore, they are well-suited for hydrostratigraphy extrapolation (Francés et al., 2014), although not for quantitative hydrogeological parameterisation. The MRS also provides hydrostratigraphic assessment (Lubczynski and Roy, 2003; Muller-Petke et al., 2011), even more reliable than other hydrogeophysical methods, but its use for hydrostratigraphy assessment is not time-wise optimal, considering the low number of survey locations that can be done per day, particularly in low signal-to-noise ratio environments (Lubczynski and Roy, 2004). To our knowledge, gravity measurements have neither been used for hydrostratigraphy assessment nor for hydrogeological system parameterisation.

## 5.3. Boundary conditions

Numerical models involve external and internal boundary conditions. Very often, those models are bound from outside by no-flow boundaries assigned along catchment water divides delineated by RS-defined DEM. The lithological boundaries, such as contacts between permeable and impermeable rocks as well as fault lines, can be delineated by multi-spectral RS, airborne or drone-borne imaging, if not entirely hidden in the subsurface. If, however, these structures are hidden, efficient hydrogeophysical methods such as EM in conductive and GPR in resistive environments are recommended.

Internal boundary conditions of models are mainly represented by surface water bodies such as rivers, streams and lakes. Temporally variable estimates of stages of surface water bodies presented in section 3.4.1 and discussed in 5.4, can serve either as time-varying specified-head boundary conditions or as observations to constrain calibrated state variables at the boundary. The temporal variability of stages can also lead to river flow and lake volume estimates. If river cross-sections and flow rating curves are known, river discharges can be deduced from
and must not be disregarded in IHMs, particularly as there are nowadays RS and hydrogeophysical tools to estimate it even in environments without in-situ measurements. Such estimates are by far more reliable than arbitrary interception assumptions, especially those neglecting interception.

Field measurements of interception (see 3.1.2) are tedious and require lots of gauges, especially in areas with complex vegetation patterns and composition; as such, they are accurate only at the plot scale. At the catchment (or coarser) scale, if interception measurements per each plant species are available, spatiotemporal interception can be scaled up by applying very high-resolution RS (Hassan et al., 2017), but also by airborne or drone-borne hydrogeophysics. However, if field interception measurements are not available or available only for some species, a use of RS (Cui and Jia, 2014; Cui et al., 2015; Gebremedhin et al., 2023) or drone-borne (Bolanos-Sánchez et al., 2021) driven analytical interception model, such as the Gash model, is recommended.

Most IHMs require as input spatiotemporally variable PET forcing (see 3.1.3). Although the PET definition is quite well-accepted and used, there is no consensus on the method for its estimation. For the PET, after Dassargues (2018), we recommend PET = ETo * Kc following the FAO approach because: i) the FAO ETo is standardised (Allen et al., 1998) and can be well defined not only by ground measurements but also by RS estimates (Gebremedhin et al., 2022); ii) the Kc for different crops is already well defined (Allen et al., 1998; Pereira et al., 2021a; Pereira et al., 2021b); iii) the Kc of natural vegetation can be estimated from spatiotemporally variable indices, such as e.g., from the NDVI, which nowadays can be derived not only from RS and airborne hydrogeophysics but also at very high resolution using Unmanned Aerial Systems (Jiang et al., 2020).

## 5.2. Parameters

IHMs require parameters of all hydrological domains. Satellite RS provides land surface and shallow soil characterisation but is of little or no use for subsurface hydrostratigraphy and parameterisation. Optical RS imagery provides useful information for setting up an IHM, such as catchment boundaries, slopes, rivers and water bodies' locations, etc. while multi-temporal optical data allows monitoring of important land surface changes (land use-land cover change, temporal variability of vegetation indices, etc.). Surface topography (DEM) is a very important IHM contribution of RS (section 3.2), drone-borne (section 4.2) and airborne (section 4.3) hydrogeophysics. However, one aspect of DEM is often underrated in setting and calibrating hydrological models: elevation accuracy. The accuracy of elevation greatly influences the reliability of IHMs because the estimates of hydraulic heads and stages of surface water bodies, directly depend on the precision of the elevation estimates. For example, how precise is a hydraulic head calculated from field measurements of water table depth with +/- 1 cm accuracy and RS elevation estimate with +/- ~10 m? Hence, standard RS DEM products with elevation accuracy of several meters must not be used for referencing hydraulic heads and stages of surface water bodies. The best elevation accuracy can obviously be obtained on the ground with geodetic, differential GPS providing accuracy of +/-1 cm, or by cumbersome, but also precise, classic geodmetrical equipment. A good option is also airborne or drone-borne green LiDAR or radar altimetry (see 4.2 and 4.3), especially in areas not accessible for ground-based elevation measurements. Considering satellite data, to our knowledge, the only comparable elevation accuracy is provided by the data of the ICESat-2 satellite mission, but its data availability declines from the poles towards the equator (see 3.3).

Parameterisation of shallow soils refers to soil type and texture, which often leads to the development of pedo-transfer functions. The microwave ability of soil penetration has been widely used for that purpose. Currently, the microwave Sentinel-1 operating in C-band (5.4 GHz) is the most suitable for that purpose, but in 2024, very high spatial resolution NISAR (NASA-ISRO Synthetic Aperture Radar) will be

launched, which will acquire data in S-band (3.2 GHz) and L-band (1.25 GHz). Soil type and texture can also be adequately provided by airborne and drone-borne hydrogeophysical platforms (section 4.2 and 4.3). In that respect, excellent results have been obtained from multi-temporal airborne microwave data acquired at 1.3 GHz and in different polarisations by Marzahn and Meyer (2020) as validated at the TERENO test-site "North-Eastern German Lowland Observatory", resulting in a mean RMSE of only 2.42 (Mass-%). The main disadvantage of microwave soil parameterisation is the limited signal penetration, which even for L-band is in order of 3–5 cm (Xu et al., 2014).

Parameterisation of the medium and deep subsurface is optimally done by combining hydrogeophysical and borehole data, while RS contribution is marginal. The hydraulic properties of the unsaturated zone can be investigated by surface hydrogeophysics (Francés and Lubczynski, 2011; Gallistl et al., 2022), but also more efficiently by airborne and recently also by drone-borne EM and GPR methods, which can also be used for spatial extrapolation of borehole data and surface hydrogeophysical surveys. Considering quantitative hydraulic parameterisation of the subsurface (unsaturated and saturated zone), the surface hydrogeophysical method of MRS is the only method providing quantitative estimates of hydrogeological parameters such as aquifer transmissivity, hydraulic conductivity and aquifer storage coefficient.

Hydrostratigraphy is an important element of conceptualised 3D pictures of the subsurface (Enemark et al., 2019). As such, it is indispensable in all constructed IHMs. 3D hydrostratigraphic models involving hydrostratigraphy and parameterisation, can be built based on borehole logs if they are sufficiently available. However, in many projects, the density of borehole logs is insufficient. Hydrogeophysical methods can provide supplementary hydrostratigraphic information to fill data gaps between the boreholes, although with different contributions, accuracy, and efficiency in different media (Enemark et al., 2020). The EM, ER, GPR and seismic reflection methods are affected by non-uniqueness mainly because the relationships between their measured physical signal and hydrogeological interpretation are proxies, confusing different hydrogeological media of similar physical properties, but they are operationally efficient. Therefore, they are well-suited for hydrostratigraphy extrapolation (Francés et al., 2014), although not for quantitative hydrogeological parameterisation. The MRS also provides hydrostratigraphic assessment (Lubczynski and Roy, 2003; Muller-Petke et al., 2011), even more reliable than other hydrogeophysical methods, but its use for hydrostratigraphy assessment is not time-wise optimal, considering the low number of survey locations that can be done per day, particularly in low signal-to-noise ratio environments (Lubczynski and Roy, 2004). To our knowledge, gravity measurements have neither been used for hydrostratigraphy assessment nor for hydrogeological system parameterisation.

## 5.3. Boundary conditions

Numerical models involve external and internal boundary conditions. Very often, those models are bound from outside by no-flow boundaries assigned along catchment water divides delineated by RS-defined DEM. The lithological boundaries, such as contacts between permeable and impermeable rocks as well as fault lines, can be delineated by multi-spectral RS, airborne or drone-borne imaging, if not entirely hidden in the subsurface. If, however, these structures are hidden, efficient hydrogeophysical methods such as EM in conductive and GPR in resistive environments are recommended.

Internal boundary conditions of models are mainly represented by surface water bodies such as rivers, streams and lakes. Temporally variable estimates of stages of surface water bodies presented in section 3.4.1 and discussed in 5.4, can serve either as time-varying specified-head boundary conditions or as observations to constrain calibrated state variables at the boundary. The temporal variability of stages can also lead to river flow and lake volume estimates. If river cross-sections and flow rating curves are known, river discharges can be deduced from