We introduce a new approach using computer vision to predict the land surface displacement from subsurface geometry images for Carbon Capture and Sequestration (CCS). CCS has been proved to be a key component for a carbon neutral society. However, scientists see there are challenges along the way including the high computational cost due to the large model scale and limitations to generalize a pre-trained model with complex physics. We tackle those challenges by training models directly from the subsurface geometry images. The goal is to understand the respons of land surface displacement due to carbon injection and utilize our trained models to inform decision making in CCS projects. We implement multiple models (CNN, ResNet, and ResNetUNet) for static mechanics problem, which is a image prediction problem. Next, we use the LSTM and transformer for transient mechanics scenario, which is a video prediction problem. It shows ResNetUNet outperforms the others thanks to its architecture in static mechanics problem, and LSTM shows comparable performance to transformer in transient problem. This report proceeds by outlining our dataset in detail followed by model descriptions in method section. Result and discussion state the key learning, observations, and conclusion with future work rounds out the paper. CCS is an important emissions reduction technology to alleviate global warming and climate change [22]. In CCS, the supercritical CO 2 is injected into the deep brine aquifers or abandoned oil reservoirs for permanent storage. However, the possible risk is the surface upheaval induced by excessive fluid pressure, which is possible to damage the structures on the ground and influence people's lives [15]. Understanding and predicting the spatial distribution and temporal change of surface deformation is important to avoid the possible hazards. However, its realization is challenged by the large spatial scale and long time scale of the project, which makes the direct numerical simulation computationally prohibitive. As the traditional numerical simulation usually requires extensive computing resources [29,1,14,9,23], the neural networks may provide a less accurate, but more efficient and sufficient solution beneficial for the model-based decision making and planning In this project, we are motivated to apply deep learning models to predict the land surface displacement under different subsurface geometries in a coupled geomechanics and flow problem. We start from the static mechanics model, where the fluid pressure is treated as the prescribed traction boundary conditions. For prediction in terms of the static mechanics model, the input to our algorithm is a geometry image. We then use different models (CNN, ResNet, ResNetUNet) to output a predicted displacement image across the domain of input geometry image. ResNetUNet is designed by us to combine the advantages of UNet and ResNet architecture to map different subsurface geometries to the distribution of the 2D displacement field. From the ResNetUNet in static mechanics model, we continue to the LSTM and transformer model to predict the temporal change of surface displacement. The input to our algorithm is a geometry image along with five frames of displacement curves, which is in format of a short video. We then use LSTM and transformer to output predicted displacement curve at the next time step. LSTM is a typical algorithm for time series predictions, and we are motivated to test transformer because of its self-attention mechanism. It provides the possibility to predict the long term dependency by contracting any two locations to a constant, so that it overcomes the issue of information loss in sequence computing in LSTM [32]. Deep learning approaches have been proved to be effective on predicting the static and time-dependent dynamics of a system under the governing physics principles as partial differential equations (PDE). The current existing deep learning approaches are of two types as physics-based neural networks (PINN) [31,17,33] and the simulation data driven approaches [25,26,30,18]. The difference is whether the direct simulation data is needed to train the model. PINN method obtains the approximate solution by minimizing the residual in physics governing equations. In contrast, common numerical discretization methods use finite difference method (FDM) and finite element method (FEM). PINN is a mesh-free method to employ automatic differentiation to handle differential operators. The existing PINN model for coupled flow and geomechanics problem [4] is limited to one dimensional problem, and it ignores the influence of material heterogeneity. Besides, it only predicts the excess pore pressure in flow side, but not the surface displacement from geomechanics side. The other type of deep learning approaches is simulation data driven method, which takes the solution from simulation directly as training data. The end-to-end image based approaches take both input and output, from figures with the model geometry as input to the contour plot showing the distribution of different physical quantities (fluid pressure, displacement) as output. The ResNet CNN was applied to predict stress field of a 2D cantilever beam under different boundary conditions [18]. Nevertheless, it is limited to static solid mechanics problem ignoring the coupling between fluid flow and mechanics response. Besides, it only deals with homogeneous model, which is not the case in reality. The R-U-Net model with LSTM is used to predict the temporal and spatial distribution of surface displacement [25], fluid pressure and saturation in the 3D model. However, it is restricted to predict one specific geological formation setting. In order to predict the temporal change of the surface displacement under different geological formation settings, a more advanced model is needed to capture the dynamics changing with time. Transformer and self-attention models [28], which were originally developed in NLP, shows its ability to outperform other methods by learning longer-term dependencies without recurrent connections. Transformer has been used to predict the representative dynamical systems such as 2D fluid dynamics and 3D reaction-diffusion dynamics [18]. The model shows its effectiveness to capture the varying behavior of different physical phenomenon evolution. Therefore, transformer is chosen to predict the transient behavior of surface displacement under different subsurface geometry setting in our project as well. In the section, we demonstrate the dataset used in the project.
The dataset includes two parts as the static model and the transient model. The input dataset are 2D geometries of geological layers , which consists of two heterogeneous materials as the hard rock with high permeability and soft shale with low permeability; and the label is surface displacement contour under specific loading. The geometries are sampled according to the real underground stratum. The size of the domain (length and width of the rectangle), the thickness, and the depth of shale layer is fixed, while the variable is the dip angle of shale layer. The data of surface displacement field is generated from simulation results in FEniCS (an open-source computational platform using Finite element Method (FEM)) [2]. The model setting used in the problem is depicted in The geometries are generated as the input of machine learning model by PIL.draw modulus, and are depicted in the Fig. 1, labeled as input image. The material characteristics of the shale and rock are distinct, thus different layout of the material will yield different results. We generate 10000 groups of data by sampling the dip angles from 0  to 45  uniformly for the training purpose. Different dip angles correspond to specific geometries, which impact the result of displacement.We split the entire dataset with 95% for training and 5% for validation. Next, we generated 5% of new samples for testing. We consider vertical displacement u y as the label, which is defined on each vertices of the finite element mesh [10], and can be considered as the displacement at points uniformly sampled over the domain. The vertical displacement (u y ) corresponds to the upheaval of the stratum. The FEniCS program provides the following API: Layout of problem = u y across the domain (1) As for the labels (ground truth), we use a numpy 1D array to store the vertical displacement of 1250 points from a uniform 50 x 25 points sampled over the domain. The dataset is demonstrated in Tab. 1. For the mathematical interpretation, refer to the Appendix A. Based on the static model, the next problem we deal with is the coupled hydro-mechanical simulation, which takes the pressure of fluid into account to better simulate the real situation. The coupled hydro-mechanical simulation is widely applied in the industry, and is well-known for dealing with consolidation case, which is a transient (time-dependent) problems. We use the same kind of geometry as what we generated for the static case, but employ another set of equations to simulate the physics process. The mathematical interpretation can be found in Appendix A. In the static case, we calculate one step of displacement across the domain, given one input geometry. However, the transient simulation will yield multiple images at different time steps given one single geometry. In addition, we focus on the vertical displacement on the ground, which is the interest in industry. In this project, we sample 40 points on the ground surface, and capture their vertical displacement. The dataset is illustrated in Tab. 2. To further prepare the data for the time series prediction model, we conduct the min-max scaling for the data by Also, we create the data sequence for the transformer and LSTM model later by using continuous five data points to predict the next data point, like where each u y,t represents a 1D array with size of 40. We adopt sliding windows to generate the sequence, thus generating 89400 pieces of displacement sequence. Similarly, we split the dataset with 95% for training, 5% for validation, and generate new samples (5%) for testing. In this project, we deal with two different computational problems, the static mechanics and the transient mechanics. We developed three models, which includes a well-designed combination of ResNet and UNet model to predict the static displacement, and the LSTM and transformer model that base on the pretrained ResNetUNet model. We are motivated to start from CNN models by literature [18], and further test ResNet and UNet from a similar deep learning problem to handle heterogeneity in input images [8]. The CNN and ResNet network reduces the dimension of images layer by layer. It loses information while collecting important features from input. Meanwhile, the deep neural network constructed by the ResNet can be utilized to decode the physical Eq. 6, which is a complicated procedure. Additionally, the heterogeneous distribution of material impacts the final result considerably. Inspired by the segmentation in computer vision [19], we use UNet to capture the heterogeneous distribution of the pixel/material. As the baseline model, we use the most basic CNN model, using continuous layers of convolution blocks which comprises a conv2d, BatchNorm2d and relu. Finally, we flatten the convolution to get a fully connected layer with 1250 dimensions for output. For the other baseline model, ResNet, we use the resnet18 model in Pytorch and modify the last layer to have the convolution with width, height and channel as 50 x 25 x 1, which will be flattened as the fully connected layer afterwords. ResNetUNet is the model designed by ourselves to solve the current problem. Compared to the conventional UNet [19], we modify the encoder part by using the ResNet for downsampling the image to better capture the physical equations. The upsampling is the same as the traditional UNet [34]. The architecture of the ResNetUNet is illustrated in the where m is the number of examples. Note that we use the MSE for the loss function, but use both MSE and maximum absolute error (MAE) for the evaluation metrics. In this part, we use both the LSTM and transformer to predict the evolution of vertical displacement, the interest in industry, given the specified material distribution as an image. We adopt two models, LSTM and transformer, to perform the task, which is similar to an image caption work [24,16], except that we use the displacement as the sequence rather than the words embedding. The combination of CNN and LSTM have been proved effective for the image caption work [24].
The CNN can be used to extract the features of the input image to feed into the LSTM network flow together with the sequence to be processed. In this project, we use a pretrained model in the static case to process the input geometry to get the features, which can embody the heterogeneous distribution of the material. The pretrained model used here is the ResNetUNet model trained in the static case, which has been proved effective for the feature extraction. Then, we add a linear layer for enhancing expressivity of the model, and processing the features further to enter the LSTM network. The linear layer has the architecture of in_features=1250, and out_feature=40. For LSTM model, we use five steps of data points as the input to predict the next single step. Since the dimension of the sequence is 40, the number of points sampled on the ground surface, thus, we use the input size of the sequence as 40. The LSTMCell model, which includes input modulation gate, input gate, forget gate and output gate [6], in Pytorch is adopted, and it is initialized by setting the input size as 40, and the hidden size as 4. We stack four LSTMCells to construct the model. The main architecture is depicted in The architecture of this model is inspired by the image caption model of LSTM [21]. The difference lies in that we use the continuous value of the displacement as the data, rather than the softmax classification of the word embedding. Thus, the loss function for the two tasks is also distinct. Based on the accuracy requirements from the industry, the dimensions of the sequence data in our engineering application are approximately the same as the dimensions of words in natural language processing. The similarity in the two models provides the confidence to build the machine learning model based on the image caption model in NLP. For the transformer model, we use the pretrained ResNetUNet to extract the image features as well, and feed it into a linear layer for secondary process. The output of the linear layer then enters the transformer encoder. As for the decoder input, we use five consecutive sequences as the input and the next five consecutive sequences as the output. For example, we use u y,t1 , u y,t2 , u y,t3 , u y,t4 , u y,t5 as the input to predict the output as u y,t2 , u y,t3 , u y,t4 , u y,t5 , u y,t6 . The input linear layer here has input features with size 1250, and out features with size 40, which reduces the dimensions of the output of pretrained model to 40. For the encoder part, we design another linear layer to convert the dimension from 40 to 16 to optimize the memory. The output then enters the positional encoding layer to take the order of the sequence into account. Then, we employ the nn.TransformerEncoderLayer in the Pytorch to assemble the encoder. The structure of the encoder is standardized from the attention paper [27]. The parameters we set for the encoder layer are d_model=16, nhead= 8. We stack four nn.TransformerEncoderLayer to form the encoder layer. For decoder, we use nn.TransformerDecoderLayer, which have the structure of d_model=16, nhead= 8. And we stack four of them to set up the decoder layer. At the end of the transformer layer, we add a linear mapping to map the dimension inner the decoder to the sequence dimension. The main architecture of the transformer is illustrated in We have studied a baseline model to predict the displacement images based on the static mechanics model. The baseline models (CNN and ResNet) are compared to a ResNetUNet model for discussion. Next, we implemented LSTM and transformer as two different models for transient mechanics model, because this is a problem related to videos. The deep learning model prediction is compared with physics modeling result, that is treated as the ground truth in this project. Due to the sufficient amount of training samples, we do not include cross validation in this work. We selected CNN as a baseline model because it collects important information from input images as features by filtering through all pixels. We have found the different subsurface geometry has strong correlation to the displacement images, and that motivates us to try CNN first to examine its performance. We generated a total of 10,000 samples and then selected 95% for training and 5% for validation. The mini batch size is 32 because this fits our CPU/GPU memory, and we prioritize the model convergence speed for less number of epochs. Although a smaller mini batch size could be less efficient in gradient calculation terms [5], a large batch may degrade the quality of the model because it tends to converge to sharp minima [12]. We tuned the hyperparameters including learning rate, optimizer, beta1, beta2, and weight decay to expect decrease of validation loss values versus the number of epochs. We also compare the training loss versus validation loss to avoid underfitting or overfitting problems. Therefore, the learning is tuned to be 5e-6, and the weight decay is 1e-5. We finally use Adam as our optimizer, because it requires less memory with large number of parameters. More importantly, the hyperparameters have intuitive interpretation, and it requires less tuning with low computational costs [13]. ResNet is the second model we investigated for static mechanics model prediction. The motivation for ResNet is high accuracy with moderate efficiency. It has deeper neural network structure so that it is capable to express more complex models. We set the same mini batch size as well as the number of epochs for training, and we follow a similar procedure to determine the hyperparameters, including learning rate, beta1, beta2, weight decay, and so on. The training loss and validation loss change versus epochs helps to determine the hyperparameters. The goal in this part is to develop a deep learning model to accurately predict the displacement images.
Therefore, it is important to visually see the predictions versus ground truth side by side. Figure 7 shows the comparison of performance based on different models. The color represents different values of displacement as a physical parameter according to the color bar, and the unit is in meter. From the ground truth image ), we observe the largest displacement happens close to the point where we apply the load. And the displacement change is getting smaller when getting farther away from the point with an applied stress. The reason to see this contour is not symmetric is that we consider heterogeneity in this system. This means the material of the subsurface rock is not the same at all different locations. The motivation comes from the subsurface engineers to characterize cases that could better represent the real world scenarios. Our ResNet prediction ) illustrates a general trend of the displacement image on the right side of the It is motivated to further improve the work by introducing ResNetUNet model because UNet includes an upsampling procedure that includes all information from each step of the downsampling process. This tends to avoid the information loss and provide more accurate predictions. We follow the same design of mini batch as well as the selection of optimizer to be Adam for this new model. And the same criterion is applied to determine hyperparameters including learning rate, beta1, beta2, and weight decay. We tuned the learning rate to be 2e-6 in order to be sufficient. And we kept the other hyperparameters the same as CNN and ResNet models. It is essential to evaluate the transient mechanics model, because this is closer to the real world problems. The physics tells us the displacement could change over time when applied a constant load on a system having porous media with fluid inside it. Therefore, the problem changes from a static image to a time series video problem. We are motivated to use LSTM because it handles time series predictions well. And the transformer is getting more and more influential, so we decided to compare LSTM and transformer models. LSTM model has shown promising potential in time series prediction in various field [20,7]. However, the current data prediction mainly refers to the low-dimensional data, which is much lower than the case we deal with here, with dimension of 40. For the training of the LSTM model, we tune parameters including the number of LSTM cells, batch_size, learning_rate. Similarly, we use an Adam optimizer for the gradient descent. To avoid the overfitting and underfitting, we finally stack 4 LTSM cells, and we choose the batch_size as 1024 to facilitate the execution in GPU. The learning_rate is chosen as 0.001. The validation loss history is shown in This paper demonstrates a state-of-the-art approach to leverage computer vision intelligence for geomechanics computations. We generalize our problem for any input geometry images, and our trained models are no longer restricted to a specific case or field. We observe the ResNetUNet model stands out compared to CNN and ResNet models for static displacement image predictions. The particular architecture of ResNetUNet ensures its excellent performance by passing the information from encoding to the decoding process. In comparison, CNN and ResNet models lose part of the information, although they successfully extract key features from input. The LSTM and transformer models work well for our transient part of work. We observe training a model for video problems requires larger RAM space due to more parameters involved. For all the models, we used both MSE and MAE as the way to quantify loss, and the comparison of training and validation loss indicates no underfitting or overfitting concerns. LSTM model shows less computational costs during the training, but we need to sacrifice prediction accuracy slightly. Both LSTM and transformer models indicate strong potential to solve the complex coupled solid mechanics and fluid flow problems, which is critical for carbon storage as one of the key components to mitigate global warming and protect our climate. For future work, we will continue the path and generate more complicated input geometry images by using Generative Adversarial Networks (GAN) to test the performance of our models. The goal is to improve our model to be robust for land surface displacement predictions to inform decision making for carbon storage projects. In addition, we would like to launch our alpha version of web application for users to access our model, so that our work will contribute value to our society. Get access to our source code in Github (Click Here). For time series prediction, refer to the demo (Click here) to view the evolution of the displacement. where  q and  p are the essential and natural boundary for the pore pressure, and the stress  is expressed as  =    pI,   = C e : (u + u). In addition to the symbols used in the static part, we have the pore pressure of fluid defined over the domain as p, and the fluid flux q = /p, which is Darcy's law that associate the flux with the pore pressure by permeability  and dynamic viscosity . The authors thank Prof. Anthony Kovscek and Prof. Ronaldo Borja for valuable suggestions and advice in terms of carbon storage and rock mechanics. We acknowledge the insightful discussions with Ruohan Gao to improve our work on both image and video problems, and we also thank the support from Manasi, Sumith to formulate our work. The authors thank Prof. Anthony Kovscek and Prof. Ronaldo Borja for valuable suggestions and advice in terms of carbon storage and rock mechanics. We acknowledge the insightful discussions with Ruohan Gao to improve our work on both image and video problems, and we also thank the support from Manasi, Sumith to formulate our work.
For the static model, the mathematical formulation is: Formulating the problem set here further in a mathematical way, we solve a partial differential problem stated as: In terms of the domian , the boundary of the domain is denoted as  =  u   t , where  u   t = , and the partial differential equations we are solving can be expressed as: where u refers to displacement of the domain, g is the gravity acceleration, and the stress  above is expressed as  = C e : (u + u), where C e is the stiffness matrix. For the transient model, the mathematical formulation is: 8a)   u +   q = 0, p = p on  p , q = q on  q , (8b) For the static model, the mathematical formulation is: Formulating the problem set here further in a mathematical way, we solve a partial differential problem stated as: In terms of the domian , the boundary of the domain is denoted as  =  u   t , where  u   t = , and the partial differential equations we are solving can be expressed as: where u refers to displacement of the domain, g is the gravity acceleration, and the stress  above is expressed as  = C e : (u + u), where C e is the stiffness matrix. For the transient model, the mathematical formulation is: 8a)   u +   q = 0, p = p on  p , q = q on  q , (8b)