## 2 Materials and Methods
### Workflow
The number of the selected pixels is far less than the full dataset size (less than 0.2%), but sufficient to take into account the main trends characterizing the input signals. The remaining larger portion of pixels is used to test the BNN. The BNN updates the weight and bias values according to Levenberg-Marquardt optimization minimizing a combination of squared errors and weights. Then, it determines the correct combination so as to produce a network that generalizes well. The BNN architecture consists of three layers. The standard structure that is used for neural network in pattern recognition applications is a two-layer feedforward network, with a sigmoid transfer function in the hidden layer, and a softmax transfer function in the output layer [78]. In particular, the softmax function is a normalized exponential function taking N-dimensional input and producing an N-dimensional output in the range (0, 1) whose values follow a probability distribution, i.e., their sum is 1. The sigmoid function is a monotonic S-shaped function taking an N-dimensional input and squashing the N-dimensional output in the range (0, 1) and their sum along classes is not 1 [78]. The first layer (input) has as many nodes as the input parameters (in our case, 13). The last layer (output) has two nodes, indicating the classes \"recent lava flow field\" (foreground) and \"outside of the recent lava flow field\" (background), with possible values 0 (not belonging to the class) and 1 (belonging to the class); the output layer nodes use a softmax activation function. A single hidden layer is interposed between the input and output layers, composed of 10 neurons with a sigmoid activation function. During the training phase, three pixels for each cluster are manually labelled to guide the following BNN training. We considered three pixels per class because they are enough to get satisfactory performance while minimizing a priori-knowledge, i.e., number of pixels to label. In particular, the three pixels closer to the medoid of each cluster are automatically selected for labelling; however, if cloudy/snowy/plume pixels are present, the expert should override the automatic choice, and label at least one of these pixels for each affected cluster, in order to avoid misclassifications. Thus, the BNN training dataset is made by three pixels for each of the \(n\) clusters, resulting in a \(3n\times 13\) BNN input matrix. The target BNN array is \(3n\times 2\), where each target pixel has value 0 or 1 for each of the background (outside of the lava flow) and foreground (inside the lava flow) classes, depending on the binary choice from the expert. Since the number of training pixels belonging to the lava flow is generally smaller than the number of training pixels outside of the flow, the user only has to identify the foreground pixels and the others are automatically labeled as background pixels. Applied to an image, the BNN will produce as output an (N \(\times\) M) \(\times\) 2 array with the classification of each pixel. From this, we generate an N \(\times\) M matrix with value 0 for pixels that are classified as background, and 1 for pixels that are classified as foreground (belonging to the \"recent lava flow field\" class). A final post-processing stage removes small isolated groups of pixels from the foreground of the 2D binary map. We assume that during an eruption, new lava will occupy a contiguous area within the flow field. In image processing, this operation is called image opening and it allows filtering out morphological noise, i.e., connected components having less than a certain number of pixels. Here, the connectivity is set to 4 and the object dimension to remove is set to less than 10 pixels. The workflow leading from the input data to the mapping of cooling lava flow fields is summarized in Figure 3. The entire code has been implemented in the MATLAB 2019a environment.