Physics-informed machine learning for understanding rock moisture dynamics in a sandstone cave

Abstract


Introduction
Water movement in the unsaturated zone is a fundamental component of the hydrologic cycle regulating the atmosphere, the hydrosphere and the lithosphere (Arora et al., 2019;Brubaker and Entekhabi, 1996;Lu and Likos, 35 2004;Tindall et al., 1998).Although there are abundant studies on water movement in various scales of unsaturated soils (Larson et al., 2022;Schoups et al., 2005;Vereecken et al., 2014;Vinnikov et al., 1996;Yu et al., 2016), much less attention has been paid to water in unsaturated rocks.In a recent study, Rempe and Dietrich (2018) defined water stored in unsaturated rocks as rock moisture and pointed out that rock moisture is a hidden component of the terrestrial hydrologic cycle critical to ecosystems and weathering processes.Due to the ubiquitous occurrence of 40 precipitation infiltration through unsaturated rocks, infiltrating precipitation was found to be the main source of rock moisture (Rempe et al., 2018;Sass, 2005).In fact, as early as in the fourth century BC, hypothesized that atmospheric water vapor could penetrate into rocks in caves with low temperature and condense into liquid water (after Meinzer, 1934).Due to occurrence of hidden water in the form of rock moisture, many stone heritages inside caves have suffered from weathering (Auler and Smart, 2004;Camuffo, 1998;de Freitas et al., https://doi.org/10.5194/hess-2022-403Preprint.Discussion started: 5 December 2022 c Author(s) 2022.CC BY 4.0 License.
for the first time and identify the atmospheric conditions controlling rock moisture in caves not exposed to rain.
Establishing the cause-and-effect relationship between rock moisture and various atmospheric conditions is a feasible approach to identify the source of rock moisture in caves and reveal mechanisms controlling rock moisture fluctuations.Machine learning has the ability to acquire knowledge and establish the complicated 60 nonlinear relationship between variables in a vast domain (Chen et al., 2019a;Jumin et al., 2020).Although machine learning models have the ability of high accuracy prediction, they are notorious for being a black-box model.Lundberg and Lee (2017) proposed the SHAP (SHapley Additive exPlanations) values as a unified measure of feature importance, which led to a combination of accuracy and interpretability of predictions by machine learning models.In almost all applications of machine learning in the field of hydrology, the directly measured 65 meteorological factors like precipitation, temperature, radiation, humidity and wind speed are used as input variables (e.g., Barzegar et al., 2017;Fang et al., 2017;Gao et al., 2020;Lees et al., 2021;Liu et al., 2022;Xiang et al., 2019;Zhao et al., 2022).In fact, using prior knowledge stemming from physical or mathematical understanding as model inputs could improve the performance of a machine learning algorithm, which is called physics-informed learning (Karniadakis et al., 2021).

70
In this study, the classic Long Short-Term Memory (LSTM) network, which is deep learning model, is combined with the SHAP values to predict rock moisture and evaluate the relative importance of four directly measured variables (precipitation, relative humidity, air temperature and wall temperature).After excluding the possible control by precipitation, based on the physics controlling vapor condensation, two new variables derived from relative humidity, air temperature and wall temperature are used as inputs of the LSTM network, which not 75 only leads to improved prediction performance, but also improve understanding of source of water in caves.Most statues in the Yungang Grottoes were carved in sandstone caves (Fig. 1a).In the summer, water droplets with planar distribution can be occasionally observed on the walls of some caves (Fig. 1b).Although no water droplets occur in other sandstone walls, by absorbing water, the high rock moisture leads to slight changes in the color of some walls.These different forms of water are responsible for weathering of the statues, however, the sources of these forms of water remains controversial.Previous studies suggested that the possible sources of water 90 in caves include infiltrating precipitation through the overlying thick unsaturated zone (Wang et al., 2012) and condensation of water vapor onto walls (Cao et al., 2005;Huang, 2008;Huang, 2010).

Monitoring of rock moisture and atmospheric conditions
To monitor variations of rock moisture in the shallow part of a cave wall, a FDR-based ECH2O EC-5 sensor (produced by DECAGON, USA) was installed at 3-8 cm inside the north wall of Cave #9 (Fig. 2).Due to the difficulty of calibrating the actual water content in the field, the apparent rock water content measured by the sensor is directly used in the current study.As reported in previous experimental studies (Mollo and Greco, 2011;Sakaki 100 and Rajaram, 2005), there is a good linear relationship between actual rock water content and rock moisture transformed from dielectric constant.Because the purpose of the current study to identify the source of rock moisture addition by establishing the relationship between rock moisture and atmospheric conditions, the relative changes of rock water content is more important than the actual value of rock water content.Therefore, direct use of apparent rock water content does not influence the results of our study, and the apparent rock water content is 105 directly called rock water content.
Air temperature (Ta) and air relative humidity (RH) are simultaneously monitored near the monitoring site of rock moisture (Fig. 2).The wall temperature is also monitored to analyze whether the wall meets the condition for condensation.Moreover, hourly precipitation is available from a meteorological station outside the cave.

Long short-term memory (LSTM) network
The LSTM network is an improved variant of the conventional Recurrent Neural Network (RNN) and has 115 the same fundamental framework as the conventional RNN.Therefore, we first give a brief introduction to the conventional RNN.As shown in Fig. 3, in the RNN, there is an input layer, a hidden layer and an output layer.xt is the input vector at time step t, ht is the hidden state at time step t determined by both the input vector xt at time step t and the hidden state (ht-1) at time step t-1 (Zhao et al., 2017), and opt is the output of the RNN at time step t.
Mathematically, the relationship between the three layers can be written as:  When the time sequence is too long, the RNN is not sufficiently efficient to learn because of vanishing 130 gradient and exploding gradient problems (Bengio, 1994;Hochreiter and Schmidhuber, 1997).To solve such problems, Hochreiter and Schmidhuber (1997) proposed an improved variant of RNN whose hidden layer can capture the correlation within time series in both short and long term.This improved variant was named the LSTM from the cell state ct should be passed to the next step (Gao et al., 2020;Fischer and Krauss, 2018;Lipton et al., 2015).Note that in the input gate, the cell state is determined simultaneously by the the sigmoid activation function which determines the weights of the information and the tanh activation function which eliminates the bias of the 140 network, the latter of which is termed abstract cell state.(forget gate)   Although purely data-driven models may have a high accuracy of prediction, most machine learning models cannot extract interpretable information and knowledge from this data deluge (Karniadakis et al., 2021).The performance of machine learning models could be improved through providing physical prior knowledge.In this 160 study, we simply integrate physics controlling vapor condensation into the input variables to improve the performance of the LSTM model.

Model interpretation and evaluation
To interpret the performance of a machine model, Lundberg and Lee (2017) proposed the SHAP (SHapley Additive exPlanations) explanation method, which is based on game theory (Štrumbelj and Kononenko, 2014). 165 The Shapely value of every input variable represents the contribution of it on the prediction, and the importance of each input variable is clarified by comparing model performances with and without it.The formula for calculation of the Shapely value is where ϕi is the contribution of variable i; F is the set of all input variables; v(S∪{i}) is the result of a model trained 170 with the variable i, and v(S) is the result without the variable i, so the difference between them represents the effect of feature i on the model prediction.This method requires retraining the model on all feature subsets S ⊆F (Shapley, 1953).
To assess the accuracy of prediction by the LSTM network, we use the statistical metrics of Nash-Sutcliffe efficiency coefficient (NSE), mean absolute error (MAE) and root mean squared error (RMSE), all of which are 175 widely used in the literature.NSE is the ratio of the sum of the squares of the regression to the total sum of the squares, which reflects the linear fit between the predictions and observations.The closer the value is to 1, the better the linear fit.The expression of NSE is where N is the number of data, and yPred, y, and ͞ y are the predicted, observed, and mean observed value, respectively.In the north wall of Cave #9, although there is no obvious occurrence of liquid water throughout the year, there is a clear trend of seasonal variation in rock water content (Fig. 5).From April to May, the rock water content is relatively stable and maintains at around 0.013 cm 3 / cm 3 .From June to September, which correspond to the rainy 190 season with high relative humidity and high air temperature, there is a cycle of significant addition and depletion of rock moisture.From October to December, there is a trend of gradual decrease in rock water content.The cycles of precipitation, relative humidity, air temperature and wall temperature from spring to early winter have quite similar trends as the cycle of rock water content, indicating that they are possible environmental conditions leading to the fluctuating rock water content.In the summer of 2021, the rock water content has a sharp increase since 9 July, reaches a maximum value equaling 0.029 cm 3 / cm 3 on 17 July, and is maintained at relatively high values until 28 July.The high rock water content indicates that there are atmospheric conditions responsible for water infiltration or condensation of water vapor.Note that this period corresponds to the period with occurrence of water droplets in Cave #5 as shown in 205 Fig. 1b1-b3.Although there is no water droplet in Cave #9, the color of the sandstone changes slightly, indicating that a slight change in the color of sandstone is a result of the fluctuating rock water content.Therefore, we believe that the rock water content measured by FDR reflects the actual change of water content in the rock.
From mid December in 2020 to the end of January in 2021, and from late December in 2021 to the end of February in 2022, there are also significant fluctuations of rock water content (Fig. 5).This pattern of fluctuating 210 rock water content is a direct of freezing-thawing, which can be confirmed by the negative wall temperature.At the beginning of the freezing-thawing cycle, there is a trend of increasing rock water content around the sensor due to freezing-induced liquid water migration towards the wall surface with the lowest temperature.By the end of the freezing period, the rock water content reaches a minimum value of the year because most liquid water has been transformed into ice.In the two years, the minimum liquid water content is 0.009 cm 3 / cm 3 (on 16 January 2021) 215 and 0.010 cm 3 / cm 3 (on 25 February 2022), respectively.In the thawing stage, there is a trend of increasing liquid water content.
The pattern of freezing-thawing-induced rock water content fluctuations is similar to that of freezingthawing-induced soil water content fluctuations (Deprez et al., 2020;Matsuoka and Murton, 2008;Sun and Scherer, 2010;Xie et al., 2021;Yu et al., 2018), demonstrating that the FDR technique is very sensitive to liquid water 220 content in sandstone and is suitable to measure rock moisture.The fluctuating rock water content during the freezing-thawing cycle also has implications for understanding weathering processes.The decreasing low liquid water content induced by freezing indicates that rock moisture in spring and autumn belongs to movable bound water that could be responsible for chemical weathering.Moreover, the freezing of accumulated liquid water near the wall surface might cause physical weathering.Deep learning models require a large amount of data for training, as well as data sets with a longtime span 235 to ensure the mastery of complete data features.Because the period from 1 June to 1 October has the most significant trends of rock moisture addition and depletion, the hourly data during this period in the year 2020 are used to construct the training set, whereas the hourly data in the year 2021 are used to construct the test set.

The predicted results based on directly monitored variables
By using relative humidity, air temperature, precipitation and wall temperature as model input variables, 240 there is a fairly good match between the predicted and measured rock water content, with similar patterns of fluctuations (Fig. 6a).Although there is obvious underestimation of rock water content in mid and late July, and slight underestimation or overestimation in other months, the NSE is as high as 0.958, indicating that the fluctuating relative humidity, air temperature, precipitation and wall temperature can capture the major patterns of fluctuating rock water content.

The prediction results based on variables controlling vapor condensation 260
Based on the SHAP values of scheme #1, precipitation can be excluded as an input variable for the LSTM network.Among the three directly monitored variables that have contributions to rock water content, air relative humidity and air temperature determines the vapor concentration and the dew point temperature (Nguyen et al., 2014), and whether the wall temperature is below the dew point temperature determines whether vapor condensation could occur (Ferná ndez-Corté s et al., 2006;Gabrovšek et al., 2010;Li et al., 2021).Because water 265 vapor is the direct source of condensation water and whether the wall temperature is below the dew point Vapor concentration and dew point temperature are both functions of actual vapor pressure, which is determined by air temperature, saturated vapor pressure and relative humidity.For air with a temperature of T (K), 270 the formula for calculating saturated vapor pressure and actual vapor pressure are (Lawrence, 2002;Lu, 2004): v,sat 273.2 0.611exp 17.27 36 where uv,sat is the saturation vapor pressure (kPa), RH is the relative humidity of air (%), and uv is the actual vapor pressure (kPa).After obtaining uv, the vapor concentration, Cv (g/m 3 ), and the dew point temperature, Td (K), 275 can be calculated as (Lu, 2004)  As indicated in Equation 16, a higher water vapor content in the air, uv, corresponds to a higher dew point temperature, thus a higher possibility of condensation at the wall.Fig. 8 shows that the patterns of fluctuating rock 280 water content, vapor concentration and difference between dew point temperature and wall temperature (denoted as Td-Tw hereafter) in the whole non-freezing period are quite similar.Moreover, we find the period with a positive Td-Tw has a good correspondence with the period with a high level of rock moisture.In scheme #2, by using the two new variables as input of the LSTM model, the mean absolute SHAP values of Td-Tw and vapor concentration are 0.0217 and 0.0100, respectively (Fig. 7b), indicating that both variables have significant contributions to rock moisture.Moreover, the NSE of predicted rock moisture is further increased to 290 0.978 (Fig. 6b).Although the prediction still underestimates rock water content from mid July to the end of July, the MAE reduced from 0.245 in scheme #1 to 0.186 in scheme #2, and RMSE reduced from 0.01416 in scheme #1 to 0.01050 in scheme #2.In the other two time durations, the MAE and RMSE of scheme #2 also decrease obviously.
Therefore, scheme #2 has much better performance of prediction, showing that using physics-informed variables would improve accuracy of prediction. 295

The mechanism of water vapor condensation
As we illustrated in 4.2, precipitation is not directly responsible for rock moisture fluctuations, but other atmospheric conditions controlling vapor concentration and the condition of vapor condensation are directly responsible for rock moisture fluctuations in the cave.In fact, vapor concentration fluctuations are more or less related to precipitation events.As shown in Fig. 9, the vapor concentration usually begins to rise prior to the 300 occurrence of a precipitation event, and declines under the control of solar radiation after a precipitation event.
Under the control of convection and diffusion, the increased water vapor in the air could be transported into porous media.When the sandstone is dry, water vapor can be adsorbed onto the surface of the rock pores, forming an adsorbed layer as thin water film; as curved menisci begin to form under increasing relative humidity, capillary condensation occurs in the rock pores (Broekhoff, 1969;Lu and Likos, 2004;Xu et al., 1998).Both adsorption and 305 capillary condensation would lead to rock moisture addition.As shown in Fig. 9, in the summer of 2021, there are 10 stages with obvious rock moisture additions.In the majority of the 10 stages, there are lagged responses of rock moisture additions to rising vapor concentration in the air, probably due to time required for vapor convection and diffusion.
Among the 10 stages, the magnitude of rock moisture addition is controlled by Td-Tw.In stages IV, V, VII, 310 VIII, IX and X, because the dew point temperature seldom exceeds the wall temperature, the magnitudes of rock moisture additions are relatively small.In stages I, II, III and VI, there are long durations with dew point temperature being higher than the wall temperature, causing large magnitudes of rock moisture addition.However, at the beginning of these four stages, even if dew point temperature is still lower than the wall temperature, increasing vapor concentration has resulted in rock moisture addition.Therefore, although a negative Td-Tw does 315 not exclude the possible occurrence of capillary condensation, a positive Td-Tw does promote capillary

Conclusion
The source of water in the sandstones caves in the Yungang Grottoes responsible for weathering was a longstanding unresolved scientific question.In this study, we use the FDR to monitor the rock moisture in a cave, which 330 shows clear rock moisture addition-depletions cycles due to various controlling mechanisms.By using relative humidity, air temperature, precipitation and wall temperature as the input variables of the LSTM network, the predicted rock moisture well reproduced the pattern of monitored rock moisture fluctuations.Moreover, we find that precipitation has no contribution, but all other three variables have contribution to the fluctuating rock moisture.
https://doi.org/10.5194/hess-2022-403Preprint.Discussion started: 5 December 2022 c Author(s) 2022.CC BY 4.0 License.accompanying a precipitation event leads to transport of water vapor into rock pores, which is subsequently adsorbed onto the surface of rock pores and then condensed into liquid water.With the aid of the deep learning model, this study increases understanding of sources of water in caves, which would contribute to future strategies 30 of alleviating weathering in caves.

Figure 1 .
Figure 1.(a) Some caves and statues in the Yungang Grottoes; (b) The occurrence (b1-b3) and disappearance (b4) of water droplets in Cave #5 of the Yungang Grottoes in the summer of 2021.95 110

Figure 2 .
Figure 2. A photo showing the arrangements of sensors for rock moisture, wall temperature, air temperature and relative humidity in the north wall of Cave #9 tanh is the activation function which means the hyperbolic tangent performs nonlinear transformations of the inputs, U, W and V are the network weight matrices for input-to-hidden, hidden-to-hidden and hidden-to-output connections, respectively, bo and bh are bias vectors.In this study, the open-source framework TensorFlow (version 125 1.14.0)written in Python 3.7.6 is used to build and train the LSTM model.
https://doi.org/10.5194/hess-2022-403Preprint.Discussion started: 5 December 2022 c Author(s) 2022.CC BY 4.0 License.network, and its hidden layer was called the LSTM cell.As shown in Fig. 4, a LSTM cell has three gates, including a forget gate, an input gate and an output gate, to maintain and adjust its cell state ct (also called memory) and 135 hidden state ht.The forget gate determines what information should be discarded from the cell state, the input gate specifies what new information need be stored in the cell state, and the output gate determines what information
https://doi.org/10.5194/hess-2022-403Preprint.Discussion started: 5 December 2022 c Author(s) 2022.CC BY 4.0 License.MAE is the mean of the distance between the predicted and the observations, whereas RMSE is the square root of the mean of the square of the deviation between the predicted and the observations.The expressions of MAE and seasonal variation of rock moisture and atmospheric conditions

195Figure 5 .
Figure 5.The temporal variations of rock water content (θ), air relative humidity (RH), air temperature (Ta), precipitation (P) and wall temperature (Tw).The two periods in yellow correspond to the summer period with high temperature and high humidity, whereas the two periods in green correspond to the fluctuation of rock moisture caused by freezing-thawing.200

2254. 2
The performance of LSTM model with two different schemesIn the rainy season, as we pointed out in 4.1, precipitation, relative humidity, air temperature and wall temperature have quite similar trends of seasonal variation as apparent rock moisture.Apparently, they are all possible factors determining the fluctuating apparent rock moisture.Therefore, in 4.2.1, we first use all of them as input variables (scheme #1) of the LSTM model to predict rock water content, and use the SHAP values to evaluate https://doi.org/10.5194/hess-2022-403Preprint.Discussion started: 5 December 2022 c Author(s) 2022.CC BY 4.0 License. the contribution of each input variable.After excluding precipitation whose mean |SHAP value| equals 0, in 4.2.2, we use two new parameters (vapor concentration, dew point temperature minus wall temperature) calculated from relative humidity, air temperature and wall temperature as input variables (scheme #2) of the LSTM model to predict apparent rock moisture.

245Figure 6 .
Figure 6.The predicted and measured rock water content obtained by two schemes.(a) Scheme #1 using four directly measured variables; (b) Scheme #2 using two calculated variables controlling vapor condensation.Also shown are NSE of the whole time series, MAE and RMSE of three different stages.

Fig
Fig.7ashows the mean absolute SHAP value of each input variable, which represents the relative importance https://doi.org/10.5194/hess-2022-403Preprint.Discussion started: 5 December 2022 c Author(s) 2022.CC BY 4.0 License.temperature is the precondition of condensation, we use vapor concentration and dew point temperature minus wall temperature as two input variables, both of which are direct variables controlling vapor condensation.

Figure 8 .
Figure 8.The fluctuating apparent rock moisture (θ), vapor concentration (Cv), and the difference between dew 285 Preprint.Discussion started: 5 December 2022 c Author(s) 2022.CC BY 4.0 License.condensation.Following the 10 stages of rock moisture additions, we find rock moisture depletion is very sensitive to decreasing vapor concentration.Moreover, in stage III with very high rock water content, a slight decrease in vapor concentration results in a slight decrease in rock water content.Therefore, we believe that rock water content 320 measured by the FDR technique is sensitive enough to fluctuating vapor concentrations and can be applied in future rock moisture monitoring in other settings.
Wch are the 155 matrices for different connections in the network, bi, bo and bc are bias vectors.In this study, the open-source framework TensorFlow (version 1.14.0)written in Python 3.7.6 is used to build and train the LSTM model.
(8)where it, ft and ot are the vectors of input, forget and output gates at time step t, respectively, all of which have the https://doi.org/10.5194/hess-2022-403Preprint.Discussion started: 5 December 2022 c Author(s) 2022.CC BY 4.0 License.samesizes as ct and ht, σ is the logistic sigmoidal activation function, c t is the vector of abstract cell state at time step t, ⊙ is element wise multiplication of two vectors.Similar toRNN, Wix, Wih, Wox, Woh, Wcx and