Reviving the Low-Frequency Response of a Rupestrian Church by Means of FDTD Simulation

: Rupestrian churches are spaces obtained from excavation of soft rocks that are frequently found in many Mediterranean countries. In the present paper the church dedicated to Saints Andrew and Procopius, located close to the city of Monopoli in Apulia (Italy) is studied. On-site acoustical measures were made, obtaining a detailed description of the acoustics in the current state pointing out, thanks to a combination of analysis techniques, the presence of significant modal behavior in the low frequencies, causing reverberation time to be about 2 s, four times longer than in the other bands, as well as being strongly dependent on source and receiver position (with variations of about 1 s when source is moved outside the chancel). However, as the church is characterized by significant degradation of surfaces and large amounts of debris cover the floor, the original acoustic conditions can be expected to somewhat differ. Acoustical modelling can be very helpful in grasping the original conditions, but given the small dimensions of the space, conventional geometrical acoustic prediction methods cannot be applied to simulate the low-frequency behavior. Thus, the present paper proposes an application of finite-difference-time-domain (FDTD) computation to simulate the low-frequency behavior and analyze a possible reconstruction of the original state. Results showed that a very good agreement was obtained between predictions and measurements, both in terms of resonance frequencies and reverberation times that differed by less than 5%. Modal response strongly affected the acoustical conditions also in the hypothetical reconstruction of the original state, although the sound field proved to be more uniform than in the current state.


Introduction
Rupestrian churches belong to the largest and most widespread group of artificial cavities of anthropic origin, often labelled as rocky or troglodytic settlements, which can be found all over the world [1]. However, in the Mediterranean area, they can be found in nearly every country, with some of them showing a significantly higher number, as outlined by a census of rocky sites in the Mediterranean area [2], from which Italy, Spain, and Turkey appear as the richest in settlements. Apulia (and the neighboring area of Matera), Sicily and Tuscany in Italy, Andalusia in Spain, and Cappadocia in Turkey showed the highest number of rocky sites, most frequently located in areas where soft stone in combination with meteorological agents already created natural caves that were subsequently enlarged and shaped according to the needs of the occupants. With reference to religious buildings, the sites in Southern Italy originated from the spreading of Greek monks following the iconoclastic prosecution in the Eastern regions but also from local communities that often found a safe place far from the frequent aggressions arriving from the sea. In Spain, many sites were built by the Mozarabic population during the Moorish occupation of the South. In Cappadocia the churches can be found in settlements of cenobite monks that developed from 6th and 7th century after the introduction of the Christian cult [3], but it is not unusual to also find mosques in some places.

The Church Surveyed
The church of Saints Andrew and Procopius (Santi Andrea and Procopio) in Monopoli (Figures 1 and 2) is located in the Contrada L'Assunta, along the Via Traiana and was likely built after the city of Monopoli was destroyed in 1042 to fight the Normans. The dedication of the church confirms the strong Byzantine influence as Andrew the Apostle was the founder of the church of Constantinople, and Procopius, martyr of Caesarea, was the protector of the Byzantine armies. The church was at the center of a rock village which also served as a post station along the Via Traiana (later Francigena) made up of a large number of caves with two or more rooms, oil mills, and mills [23].
Acoustics 2023, 5 4 FOR PEER REVIEW presents the results of the on-site measurements, the results of the FDTD simulation in th current conditions, and the reconstruction of the original state; Section 4 provides a bri discussion of the results in comparison with other existing studies; finally, Section summarizes the major conclusions.

The Church Surveyed
The church of Saints Andrew and Procopius (Santi Andrea and Procopio) Monopoli (Figures 1 and 2) is located in the Contrada L'Assunta, along the Via Traian and was likely built after the city of Monopoli was destroyed in 1042 to fight the Norman The dedication of the church confirms the strong Byzantine influence as Andrew th Apostle was the founder of the church of Constantinople, and Procopius, martyr Caesarea, was the protector of the Byzantine armies. The church was at the center of rock village which also served as a post station along the Via Traiana (later Francigen made up of a large number of caves with two or more rooms, oil mills, and mills [23].  The facade, just like in a masonry temple, has three arched entrances. The church is divided in two parts, the "naos" occupied by the congregation, having a simple rectangular shape (about 5.6 m by 5.2 m) and the chancel, divided by means of a "templum" (i.e., a stone iconostasis). The chancel was subdivided into four squares, two of them (the farthest from the entrance) define the "bema" where the altar was located. The ceiling is flat and its current height varies between 1.90 m and 2.40 m, but the original height is likely to be greater as a large amount of debris covers the floor. All the walls were originally The facade, just like in a masonry temple, has three arched entrances. The church is divided in two parts, the "naos" occupied by the congregation, having a simple rectangular shape (about 5.6 m by 5.2 m) and the chancel, divided by means of a "templum" (i.e., a stone iconostasis). The chancel was subdivided into four squares, two of them (the farthest from the entrance) define the "bema" where the altar was located. The ceiling is flat and its current height varies between 1.90 m and 2.40 m, but the original height is likely to be greater as a large amount of debris covers the floor. All the walls were originally covered by frescos but only a few fragments may be admired today.
Among the frescoes on the walls that are still visible, there are the apostles Peter and Paul, the saints Cosma, Damiano, Eligio (with the symbols of his patronage over the blacksmiths), Giorgio, Leonardo, a Virgin in throne and a scene of the Annunciation, datable between the 13th and 14th centuries. The painted scenes are all found in the bema and transept: Annunciation, Deesis, Trinity, Crucifixion, while the saints instead occupy the walls of the naos. The wall decoration is not contemporaneous with the excavation of the crypt but was constructed two centuries later. However, it reflects the continuity of worship from the Norman age, to which Saints Peter and Paul, Leonardo and Eligius belong, for example, up to the period of the Crusades (see, for example, the fresco of Saint George with the emblem of the cross on the saddle and on the shield). Unfortunately lost is the Byzantine iconography of the eleventh century relating to the eponymous saints of the sanctuary.
From the point of view of the liturgical functions, the church had the bema closed by a first iconostasis in stone opened by two doors in correspondence with the two apses. The transverse arm of the transept separates the presbytery from the naos, a common quadrangular room, not divided by pillars, with a second iconostasis, also with two doors. A summary of the main geometrical features is given in Table 1.

Acoustic Measurement Methods
Given the location of the church in an open field, all the measurements were carried out with portable instruments powered with battery packs. An omni-directional sound source (Lookline D301) was located in two positions, one in front of the altar (Position A) and one in the congregation area (Position B), so as to simulate in this way the priest and the congregation, respectively. The source was fed by an equalized sine sweep played back by a portable music player and generated using MATLAB (v. 2021b) according to Müller and Massarani [24] so that the spectrum of the radiated sound was substantially flat from 50 Hz to 16 kHz. The duration of the sweep was kept short (about 8 s) in order to limit any potentially adverse effects due to lack of doors, which, determining significant air circulation, might compromise the linear and time-invariant hypothesis. In addition, given the conservation state of the surface finishing, longer sine sweeps might have induced excess mechanical stress, so it was preferred to keep signals short. Given the limited dimensions of the space, the signal-to-noise (S/N) ratio was sufficiently high to ensure perfectly usable impulse responses, even if a short sweep was used. Room responses were collected using a portable B-format microphone (Soundfield ST-350) connected to a multi-channel recorder (Tascam DR-680) and a pair of binaural microphones (Soundman OKM II) worn by one of the authors and connected to a second recorder (Tascam DR-07). The measurement chain was previously tested in the lab to ensure that the "open loop" settings did not create any sync problem.
All the measurements were carried out complying with ISO 3382-1 [25] standard and taking into account the guidelines for measurements in worship spaces [26], and, despite the small dimensions of the church, eight receiver positions were used to provide a detailed description of the point-by-point variations. Microphones were placed at 1.6 m from the floor, assuming that the congregation was standing during the celebration. Source and receiver locations were chosen to provide a description of acoustical conditions that could be produced by actual sound sources and heard by listeners distributed in the space. Given the small dimension of the churches, only one person stayed in the room during the measurements.
Impulse responses (IR) were calculated by deconvolving the signal used to feed the sound source and, despite a significant background noise due to birds and other natural sounds (resulting, on average, in an A-weighted sound pressure level of 45 dB), provided a minimum S/N ratio of about 55 dB over the whole spectrum of interest. The measured IRs were then processed in order to calculate the most important acoustic parameters and to investigate room resonances. In particular, in addition to monaural parameters based on the omni (W) response of the B-format microphone, lateral energy fraction was calculated using W and Y Ambisonic components (assuming X axis was aimed at the source), while inter-aural cross correlation was based on binaural responses.

Geometrical Modelling of the Space
Given the complex and irregular shape of the space, the only reasonable way to obtain a reliable geometrical model was to use 3D laser scanning. A point cloud of more than 150k elements was originally obtained using a Riegl VZ 400 scanner (with an original resolution of 5 cm), but for the purpose of the simulation, such level of detail was unnecessary and thus, after cleaning artifacts and imperfections using a specifically designed tool developed in MATLAB, the point cloud was further simplified using the open source software Meshlab through subsequent applications of the "Quadric Edge Collapse Decimation" algorithm ( Figure 3).

Geometrical Modelling of the Space
Given the complex and irregular shape of the space, the only reasonable way to obtain a reliable geometrical model was to use 3D laser scanning. A point cloud of more than 150k elements was originally obtained using a Riegl VZ 400 scanner (with an original resolution of 5 cm), but for the purpose of the simulation, such level of detail was unnecessary and thus, after cleaning artifacts and imperfections using a specifically designed tool developed in MATLAB, the point cloud was further simplified using the open source software Meshlab through subsequent applications of the "Quadric Edge Collapse Decimation" algorithm ( Figure 3).

FDTD Simulation
FDTD acoustic modelling has been applied to solve acoustic problems for a long time [27] but given the computational load it has been mostly applied to low frequencies [28], while the availability of parallel computation distributed by several GPUs fostered a gradual extension to a much wider frequency range [21,22,[28][29][30]. FDTD starts from the assumption that a generic derivation operator can be replaced by one of its finite difference forms: and hence the second-order derivative becomes:

FDTD Simulation
FDTD acoustic modelling has been applied to solve acoustic problems for a long time [27] but given the computational load it has been mostly applied to low frequencies [28], while the availability of parallel computation distributed by several GPUs fostered a gradual extension to a much wider frequency range [21,22,[28][29][30]. FDTD starts from the assumption that a generic derivation operator can be replaced by one of its finite difference forms: and hence the second-order derivative becomes: In acoustics, any propagation phenomenon is described by the wave equation expressed by: The equation is only expressed as a function of sound pressure p, but its finitedifference form would require three different time steps to be numerically solved, so in a 2D case in which T is the time step and X is the grid spacing, pressure at step n + 1 and node (x,y) will be given by: p n+1 x,y = λ 2 p n x+1,y + p n x−1,y + p n x,y+1 + p n x,y−1 + 2 1 − 2λ 2 p n x,y − p n−1 x,y where λ = c T/X, c being the speed of sound in air. However, considering that the wave equation is derived by two other fundamental equations, involving only first-order derivatives, but also including particle velocity, it is possible to find a perfectly equivalent [29,30] alternative formulation. This formulation is known as staggered Yee's grid, in which the grid of the pressure values is complemented by the grid of the particle velocity values (u). Thus, it is possible to first update particle velocity components (remembering that particle velocity is a vector quantity and, in 2D, it will have ux and uy components) and when this matrix is available, we can update the pressure at the next time step. Pressure at step n + 1 and node (x,y) will consequently be given by: x,y+0.5 − uy n+0.5 x,y−0.5 The above equations apply when non-boundary conditions are found, while at grid points close to walls, assuming a wall impedance Z, it is possible to find an update formula that, in case Z may be considered a real number independent of frequency, yields (in the 2D case) [31]: where being the absorption coefficient of the given surface. A number of alternative approaches have been proposed to account for frequency dependence [27,32], but in the present case, being the analysis limited to the low-frequency case, the proposed assumption was not considered a major limitation. In order to obtain a reliable calculation using FDTD, it is essential to properly set grid spacing X and time step T. In fact, they are both related to the maximum frequency that can be analyzed, as a minimum of five points per wavelength is usually required to prevent errors. Therefore, this means that given the minimum wavelength of interest (λ), the grid spacing should be at least λ/5 or, better, λ/10. Thus, for a grid spacing of 0.1 m, the FDTD results will be accurate up to a maximum frequency of 340 Hz. The time step is not independent of the grid spacing if the stability condition given by: (where D is the number of dimensions of the problem) is satisfied. Thus, once the grid spacing and the problem dimension are defined, time step results consequently (in the present case corresponding to a sampling frequency of 5768 Hz). Proper limitation of the highest frequency of analysis is also useful to minimize dispersion errors that result in phase velocity being different from the actual value of the medium. In the present case, limiting to frequencies below 340 Hz ensures that dispersion errors will be below 2% [29]. The previously described FDTD framework was implemented in MATLAB (Figure 4), where the geometrical model was first voxelized using a 10 cm grid spacing, and then surface properties could be assigned together with source and receiver locations. Sound sources were modeled as simple point sources and emitted signals could vary between sine waves (for modal analysis) and short pulses. As per usual, surface properties were adjusted in order to have a suitable match between measured and predicted values of the reverberation time in the present condition at frequencies of 63 Hz and 125 Hz.
It is important to point out that for a space of such dimensions (113 m 3 ), it is essential to model sound propagation by means of wave-based models because geometrical acoustics can only be effective well above Schroeder's frequency that, being equal to 2000 / , in this case was around 250 Hz. For this reason, modal response was also carefully compared to ensure that the model could realistically reproduce low-frequency propagation.

Material Characterization
In order to characterize sound absorption of the tuffaceous surface of the church, samples taken from quarries of geologically similar limestone were analyzed. Measurements of normal incidence sound absorption coefficient were carried out according to ISO 10534-2:1998 [33], using the transfer function method. As the objective of the study was to understand the low-frequency behavior of the materials, only the tube with an internal diameter of 10 cm was used, resulting in a maximum measurable frequency of 2 kHz and a low-frequency limit of 50 Hz. The emitting end consisted of an 11 cm loudspeaker sealed into a wooden case and suitably isolated from the tube structure by an elastic and protective layer. All the processing was performed by a MATLAB graphic user interface gener- As per usual, surface properties were adjusted in order to have a suitable match between measured and predicted values of the reverberation time in the present condition at frequencies of 63 Hz and 125 Hz.
It is important to point out that for a space of such dimensions (113 m 3 ), it is essential to model sound propagation by means of wave-based models because geometrical acoustics can only be effective well above Schroeder's frequency that, being equal to 2000 √ T/V, in this case was around 250 Hz. For this reason, modal response was also carefully compared to ensure that the model could realistically reproduce low-frequency propagation.

Material Characterization
In order to characterize sound absorption of the tuffaceous surface of the church, samples taken from quarries of geologically similar limestone were analyzed. Measurements of normal incidence sound absorption coefficient were carried out according to ISO 10534-2:1998 [33], using the transfer function method. As the objective of the study was to understand the low-frequency behavior of the materials, only the tube with an internal diameter of 10 cm was used, resulting in a maximum measurable frequency of 2 kHz and a low-frequency limit of 50 Hz. The emitting end consisted of an 11 cm loudspeaker sealed into a wooden case and suitably isolated from the tube structure by an elastic and protective layer. All the processing was performed by a MATLAB graphic user interface generating a linear sweep to feed the loudspeaker.
As shown in Table 2, the major difference between different types of tuffaceous stones appears above 80 Hz, where higher porosity may contribute to nearly double absorption compared to harder samples. At very low frequencies, no significant differences were observed in the measured values.

On-Site Acoustic Measurements
The analysis of the reverberation times shows very interesting differences between the medium-high frequencies, where the values settle below 0.5 s, with negligible point-bypoint variations, and the low frequencies, where the values are considerably longer and show some variability depending on the location of both source and receivers ( Figure 5). In particular, when the source was close to the altar (Source A), T30 at 63 Hz and at 125 Hz assumed the longest values (particularly at receivers 1, 2, 4, and 7), smoothly decreasing when moving closer to the entrance. When the source was in position B, T30 dropped by about 0.5 s in the same position. This fact, combined with EDT values shorter than T30 whatever the source position, suggested that no evident reverberant coupling effects between the sub-volumes may explain the observed variations. Conversely, the analysis of the time decays at 63 Hz and 125 Hz shows ( Figure 6) that decays are characterized by a staircased or "pulsating" trend, clearly evident at 63 Hz, but also appearing at 125 Hz, in particular for combinations A 01 and A 02. Such behavior is typically associated with repeated reflections (flutter echoes) or modal effects which could be better investigated by analyzing the narrow band spectra of the responses.

On-Site Acoustic Measurements
The analysis of the reverberation times shows very interesting differences between the medium-high frequencies, where the values settle below 0.5 s, with negligible pointby-point variations, and the low frequencies, where the values are considerably longer and show some variability depending on the location of both source and receivers ( Figure  5). In particular, when the source was close to the altar (Source A), T30 at 63 Hz and at 125 Hz assumed the longest values (particularly at receivers 1, 2, 4, and 7), smoothly decreasing when moving closer to the entrance. When the source was in position B, T30 dropped by about 0.5 s in the same position. This fact, combined with EDT values shorter than T30 whatever the source position, suggested that no evident reverberant coupling effects between the sub-volumes may explain the observed variations. Conversely, the analysis of the time decays at 63 Hz and 125 Hz shows ( Figure 6) that decays are characterized by a staircased or "pulsating" trend, clearly evident at 63 Hz, but also appearing at 125 Hz, in particular for combinations A_01 and A_02. Such behavior is typically associated with repeated reflections (flutter echoes) or modal effects which could be better investigated by analyzing the narrow band spectra of the responses.  Before delving into the analysis of the spectra, it is worth noticing that the short reverberation time at medium-high frequencies was likely due to the presence of openings (about 4 m 2 ) and to the particular state of deterioration of the stone which appeared very pronounced, with frescoes occupying only a small part of the surfaces and the underlying layer characterized by greater porosity and tenderness. As shown in Table 2, a soft limestone is capable of absorbing a significantly higher amount of acoustic energy, which may easily explain the observed values. In addition, the floor was completely covered by waste soil overflowed during floods, characterized, at the time of measurement, by numerous cracks and having a thickness which could be estimated to vary between 10 cm and 50 cm in some points.
In order to better understand the observed low-frequency behavior, narrow band spectra were determined for all the source-receiver combinations ( Figure 7). As expected, when source was in position A, receivers 1, 2, and 4 clearly showed a marked peak at 79 Hz, surrounded by many others appearing at 65 Hz, 73 Hz, and 87 Hz. In the other receivers, the same modes appeared but their energy content was significantly reduced. In fact, if the mid-frequency spectrum density is taken as a reference, the level variation around 79 Hz is about 30 dB between receivers 1 and 2 and the others. Before delving into the analysis of the spectra, it is worth noticing that the short reverberation time at medium-high frequencies was likely due to the presence of openings (about 4 m 2 ) and to the particular state of deterioration of the stone which appeared very pronounced, with frescoes occupying only a small part of the surfaces and the underlying layer characterized by greater porosity and tenderness. As shown in Table 2, a soft limestone is capable of absorbing a significantly higher amount of acoustic energy, which may easily explain the observed values. In addition, the floor was completely covered by waste soil overflowed during floods, characterized, at the time of measurement, by numerous cracks and having a thickness which could be estimated to vary between 10 cm and 50 cm in some points.
In order to better understand the observed low-frequency behavior, narrow band spectra were determined for all the source-receiver combinations ( Figure 7). As expected, when source was in position A, receivers 1, 2, and 4 clearly showed a marked peak at 79 Hz, surrounded by many others appearing at 65 Hz, 73 Hz, and 87 Hz. In the other receivers, the same modes appeared but their energy content was significantly reduced. In fact, if the mid-frequency spectrum density is taken as a reference, the level variation around 79 Hz is about 30 dB between receivers 1 and 2 and the others.
When the source was moved into the "naos", the acoustic energy redistributed among modes. In fact, the new absolute peak appeared around 65 Hz and 73 Hz at receivers 1, 2, and 4, while in the others the mode energy was more evenly distributed but remained at least 10 dB higher than in the frequencies above 100 Hz.
This can clearly explain what was observed in terms of reverberation time because the peak appearing around 80 Hz due to its position and magnitude influenced both the octave bands of 63 Hz and 125 Hz, causing the slower decay as observed in Figure 6 to appear in receiver 1 and 2 also in the higher octave band. In the other receivers, as the energy in the modes decreases, the reverberation time becomes gradually shorter, remaining longer in the 63 Hz band due to persistence in time of modal behavior (Figure 8a). To this purpose, modal reverberation time [34,35] was calculated for combination A-01 at the frequency of 78.2 Hz, showing (Figure 8b) that the observed value basically coincided with the octave band value. Acoustics 2023, 5 4 FOR PEER REVIEW 11 Rec 01 Rec 02 Rec 03 Rec 04 Rec 05 Rec 06 Rec 07 Rec 08 When the source was moved into the "naos", the acoustic energy redistributed among modes. In fact, the new absolute peak appeared around 65 Hz and 73 Hz at receivers 1, 2, and 4, while in the others the mode energy was more evenly distributed but remained at least 10 dB higher than in the frequencies above 100 Hz. appear in receiver 1 and 2 also in the higher octave band. In the other receivers, as the energy in the modes decreases, the reverberation time becomes gradually shorter, remaining longer in the 63 Hz band due to persistence in time of modal behavior (Figure 8a). To this purpose, modal reverberation time [34,35] was calculated for combination A-01 at the frequency of 78.2 Hz, showing (Figure 8b) that the observed value basically coincided with the octave band value. The mathematics of the axial resonance modes [20] state that resonances appear at frequency fn = (c ⋅ n)/(2L), where c is the speed of sound, L is the room dimension in the direction that is considered, and n is an integer index that defines the mode order. Thus, it can be easily found that 79 Hz corresponds to the first (n = 1) axial mode along the vertical direction assuming a height of 2.15 m (which is in good agreement with the average height of the space, although being far from a simple rectangular box), while 65 Hz corresponds well to the second (n = 2) axial mode along the two side walls (spaced by about 5.4 m), the only two large and continuous surfaces besides the floor and the ceiling.
Once the basic acoustical features have been explained by properly combining reverberation time and modal analysis, it is now possible to have a look at the other acoustical parameters (Figure 9). Given the short reverberation, very high speech intelligibility is obtained (STI = 0.82 on average) and strong frequency imbalance represented by a bass ratio (BR) equal to 2.1, although in the original conditions, with walls covered by frescos, it is likely that a more balanced condition was observed. The mathematics of the axial resonance modes [20] state that resonances appear at frequency f n = (c · n)/(2L), where c is the speed of sound, L is the room dimension in the direction that is considered, and n is an integer index that defines the mode order. Thus, it can be easily found that 79 Hz corresponds to the first (n = 1) axial mode along the vertical direction assuming a height of 2.15 m (which is in good agreement with the average height of the space, although being far from a simple rectangular box), while 65 Hz corresponds well to the second (n = 2) axial mode along the two side walls (spaced by about 5.4 m), the only two large and continuous surfaces besides the floor and the ceiling.
Once the basic acoustical features have been explained by properly combining reverberation time and modal analysis, it is now possible to have a look at the other acoustical parameters (Figure 9). Given the short reverberation, very high speech intelligibility is obtained (STI = 0.82 on average) and strong frequency imbalance represented by a bass ratio (BR) equal to 2.1, although in the original conditions, with walls covered by frescos, it is likely that a more balanced condition was observed. The perception of sound spatiality and, more generally, the binaural impression were particularly interesting. JLF and 1-IACC applied to the initial part of the impulse response showed rather high values (1-IACC = 0.67 on average) which were reduced only (as expected) in correspondence with the points placed in direct proximity to the source (receiv- The perception of sound spatiality and, more generally, the binaural impression were particularly interesting. JLF and 1-IACC applied to the initial part of the impulse response showed rather high values (1-IACC = 0.67 on average) which were reduced only (as expected) in correspondence with the points placed in direct proximity to the source (receivers 2 and 5). The large amount of surface irregularities distributed almost everywhere, together with the presence of obstacles (such as the iconostasis) all contribute to this effect, making the sound perception very "spacious" despite the small size. JLF showed extremely high values deriving from the strong contribution of lateral reflections combined with the frequent shielding of reflections coming from frontal directions.

FDTD Acoustical Simulation of Current State
FDTD simulations were carried out according to the procedures described in Section 2.4. As the finishing was relatively the same over all the surfaces, with the notable exception of the floor covered by hardened soil, and given the values measured for similar materials (Table 2), an initial value of 0.015 was used for the absorption coefficient (obtained by averaging among the 63 Hz and the 125 Hz values). As this value returned slightly longer T30 values than expected, while in acoustic modelling [22] a maximum difference of 5% (corresponding to one just noticeable difference [25]), is considered to be acceptable, step-by-step increases were applied until a value of 0.025 was reached. Under these conditions, the spatial average of T30 in one-third octave bands from 63 Hz to 80 Hz differed by only 2% from measurements, while at 100 Hz the error was 4%. At higher frequencies, the error was bigger, suggesting that a further increase in absorption coefficient was needed. By adopting an α value of 0.04, the mean error was finally reduced to 4% also in the 125 Hz band. Figure 10 shows the comparison between measured and predicted reverberation time in one-third octave bands at individual positions. It can be observed that in the lowest bands clear modal behavior appears at given receiver positions, resulting from the strong non-uniform distribution of the values. In addition, as discussed in the previous section, the strength of the first axial mode along the vertical direction clearly appears also in the simulated results, where T30 in the 80 Hz band is markedly higher than the corresponding values at 63 Hz despite the same absorption coefficient. At 100 Hz and 125 Hz the modal behavior is significantly attenuated in the simulations (although the measured values at combinations A-01 and A-02 still show some effects). Given the role played by modal response of the room, it was important to check the agreement between measured and predicted spectra. Considering that Figure 7 had already shown the combinations where a strong modal response was measured, the comparison was restricted to a subset of all the source-receiver combinations. Figure 11 shows that for combinations A-01 and A-02, the strong response around 79 Hz clearly appears, although the magnitude is attenuated compared to measurements assuming higher frequencies as a reference (which also explains why the 100 Hz T30 is shorter in the simulated   Given the role played by modal response of the room, it was important to check the agreement between measured and predicted spectra. Considering that Figure 7 had already shown the combinations where a strong modal response was measured, the comparison was restricted to a subset of all the source-receiver combinations. Figure 11 shows that for combinations A-01 and A-02, the strong response around 79 Hz clearly appears, although the magnitude is attenuated compared to measurements assuming higher frequencies as a reference (which also explains why the 100 Hz T30 is shorter in the simulated response). Small shifts in the peak frequencies can be observed, but they rarely exceed 2 Hz. Combination A-05 is modelled pretty well, with the peak around 65 Hz well visible and a slightly less defined cluster around 79 Hz but a well comparable distribution of levels. Finally, combination B-04 is considered, as in the measured spectrum a clear peak appeared at 65 Hz. The simulated spectrum presents several peaks that are mostly aligned with measurements, but the stronger peak at 65 Hz seems more attenuated (by about 15 dB), and another peak below 60 Hz appears. However, considering that measured spectra were band-limited by the inherent frequency response of the sound source and by the signal used to feed the loudspeaker, this was considered a minor problem. A further element to be considered that might explain some of the small differences observed in the T30 values and in the responses is the strong point-by-point variations that appear in consequence of the modal behavior. Taking into account the FDTD model, the sound pressure level resulting from pure tones at 65 Hz and 79 Hz was calculated for the two source positions. As shown in Figure 12, dramatic sound pressure level variations take place by simply moving the receivers a few centimeters apart. Hence, it is not unlikely that small misplacement of the receivers (both during the measurements and in the simulation setup) might explain some of the observed differences.
However, Figure 12 is quite instructive because it shows which are the positions in the church where stronger resonances appear and, consequently, where longer reverberation might be experienced. Predictably, receivers close to the side walls in the "naos" could experience strong resonances whatever the source position was, particularly at 65 Hz (which was the second axial mode along the transverse direction), but also at the other frequencies. A further element to be considered that might explain some of the small differences observed in the T30 values and in the responses is the strong point-by-point variations that appear in consequence of the modal behavior. Taking into account the FDTD model, the sound pressure level resulting from pure tones at 65 Hz and 79 Hz was calculated for the two source positions. As shown in Figure 12, dramatic sound pressure level variations take place by simply moving the receivers a few centimeters apart. Hence, it is not unlikely that small misplacement of the receivers (both during the measurements and in the simulation setup) might explain some of the observed differences.

FDTD Acoustical Reconstruction of Original State
Once the model was calibrated, it could be used to apply a few changes that could realistically correspond to the original state, allowing us to appreciate the implications for the acoustics. In the specific case, three simple changes were applied. First, the floor area was rectified assuming that the layers of debris that currently occupy the entrance and the chancel area, where the ceiling is unusually low (see dashed line in Figure 1b), could be removed. Surfaces were considered to be covered by plasters and frescos (like the limited portions still existing clearly suggest, and like it is found in many other rupestrian churches in the region), resulting in a reduced absorption coefficient of 0.015 at 63 Hz and 80 Hz, and 0.025 at 100 Hz. Finally, in order to simulate acoustics under occupied conditions, a seated audience was located along the two side walls of the "naos", according to the common practice of the time, as demonstrated by the usual presence of "subsellia" (carved benches) in other churches. This position was chosen in order to locate the absorption in the position where modal behavior might have more strongly affected the acoustic response in terms of damping. Absorption coefficients for this area were set to 0.2 [36] assuming the audience to be tightly distributed among the limited seats. All the openings were left in their actual state as no evidence of pre-existing doors or windows was found.
As shown in Figure 13, most of the main features observed in the "current" state also appear in the "reconstructed", but now reverberation is longer and the modal response at 79 Hz when the source is in A is stronger and extends well beyond the "templum". When the source is in B, the most evident variation is a drop appearing at receivers close to the entrance, possibly as a consequence of the change in room height (that is the part where the debris were thicker and after removal the ceiling height becomes about 2.7 m) and increase in the absorbing elements in the volume due to audience and increased opening surfaces (again as a consequence of debris removal). However, Figure 12 is quite instructive because it shows which are the positions in the church where stronger resonances appear and, consequently, where longer reverberation might be experienced. Predictably, receivers close to the side walls in the "naos" could experience strong resonances whatever the source position was, particularly at 65 Hz (which was the second axial mode along the transverse direction), but also at the other frequencies.

FDTD Acoustical Reconstruction of Original State
Once the model was calibrated, it could be used to apply a few changes that could realistically correspond to the original state, allowing us to appreciate the implications for the acoustics. In the specific case, three simple changes were applied. First, the floor area was rectified assuming that the layers of debris that currently occupy the entrance and the chancel area, where the ceiling is unusually low (see dashed line in Figure 1b), could be removed. Surfaces were considered to be covered by plasters and frescos (like the limited portions still existing clearly suggest, and like it is found in many other rupestrian churches in the region), resulting in a reduced absorption coefficient of 0.015 at 63 Hz and 80 Hz, and 0.025 at 100 Hz. Finally, in order to simulate acoustics under occupied conditions, a seated audience was located along the two side walls of the "naos", according to the common practice of the time, as demonstrated by the usual presence of "subsellia" (carved benches) in other churches. This position was chosen in order to locate the absorption in the position where modal behavior might have more strongly affected the acoustic response in terms of damping. Absorption coefficients for this area were set to 0.2 [36] assuming the audience to be tightly distributed among the limited seats. All the openings were left in their actual state as no evidence of pre-existing doors or windows was found.
As shown in Figure 13, most of the main features observed in the "current" state also appear in the "reconstructed", but now reverberation is longer and the modal response at 79 Hz when the source is in A is stronger and extends well beyond the "templum". When the source is in B, the most evident variation is a drop appearing at receivers close to the entrance, possibly as a consequence of the change in room height (that is the part where the debris were thicker and after removal the ceiling height becomes about 2.7 m) and increase in the absorbing elements in the volume due to audience and increased opening surfaces (again as a consequence of debris removal). Such results seem to suggest that the resulting acoustic effects might have been exploited during liturgical celebrations and singing, in particular, considering that the specific features of many Byzantine hymns include many bass notes sustained for a long time, which could well excite the resonances of the space.

Discussion
After presenting the results of the measurements and of the simulations, a brief discussion to put the results in the context of the existing literature can be developed. Among the spaces that have been surveyed by other authors, the highest similarity can be found with the rock-cut structures in Cappadocia [18]. Although the only space having a comparable volume was the Avanos dining room (114 m 3 ), and the observed reverberation times spanned over a much broader range (up to 5 s for the Hallaç Church and main hall and for Açıksaray Hall), in all of the cases a significant low-frequency imbalance was observed, with reverberation time being up to 3 times longer than at mid-frequencies. This was likely caused by the characteristics of the stone that due to its porosity caused increased absorption as frequency grew. However, no specific low-frequency measures were made, and the frequency range of analysis was limited to 125 Hz.
Similarly, among the catacombs in Southern Italy [17], it was difficult to find spaces that are entirely comparable, but the two cubicles in the "San Callisto" catacombs have several similarities in terms of room dimensions and volumes. Materials are also similar, the walls being made of tufa stone, resulting in short reverberation times usually well below 1 s. In this case, the frequency range of the measurements extends into the 63 Hz band, showing only a moderate increase in reverberation time compared to mid-frequencies, with the notable exception of the so called "double cubicle", where reverberation is longer (up to 1.5 s at 63 Hz) and characterized by point-by-point variations. Regardless, even in this study, although the "resonant" behavior in the low frequencies is supposed, no additional investigations were made.
Low-frequency resonances and modal behavior were investigated in many stone chambers and cairns in the British Isles [7]. In this case, experimental measurements were carried out using sine sweeps to find resonant frequencies and the distribution of nodes and antinodes in the space. A theoretical justification was consequently found based on the dimensions of the space, showing that in most of the cases, resonance frequencies Such results seem to suggest that the resulting acoustic effects might have been exploited during liturgical celebrations and singing, in particular, considering that the specific features of many Byzantine hymns include many bass notes sustained for a long time, which could well excite the resonances of the space.

Discussion
After presenting the results of the measurements and of the simulations, a brief discussion to put the results in the context of the existing literature can be developed. Among the spaces that have been surveyed by other authors, the highest similarity can be found with the rock-cut structures in Cappadocia [18]. Although the only space having a comparable volume was the Avanos dining room (114 m 3 ), and the observed reverberation times spanned over a much broader range (up to 5 s for the Hallaç Church and main hall and for Açıksaray Hall), in all of the cases a significant low-frequency imbalance was observed, with reverberation time being up to 3 times longer than at mid-frequencies. This was likely caused by the characteristics of the stone that due to its porosity caused increased absorption as frequency grew. However, no specific low-frequency measures were made, and the frequency range of analysis was limited to 125 Hz.
Similarly, among the catacombs in Southern Italy [17], it was difficult to find spaces that are entirely comparable, but the two cubicles in the "San Callisto" catacombs have several similarities in terms of room dimensions and volumes. Materials are also similar, the walls being made of tufa stone, resulting in short reverberation times usually well below 1 s. In this case, the frequency range of the measurements extends into the 63 Hz band, showing only a moderate increase in reverberation time compared to mid-frequencies, with the notable exception of the so called "double cubicle", where reverberation is longer (up to 1.5 s at 63 Hz) and characterized by point-by-point variations. Regardless, even in this study, although the "resonant" behavior in the low frequencies is supposed, no additional investigations were made.
Low-frequency resonances and modal behavior were investigated in many stone chambers and cairns in the British Isles [7]. In this case, experimental measurements were carried out using sine sweeps to find resonant frequencies and the distribution of nodes and antinodes in the space. A theoretical justification was consequently found based on the dimensions of the space, showing that in most of the cases, resonance frequencies varied between 95 Hz and 120 Hz, slightly above the resonances observed in the present work, likely as a consequence of the smaller dimension.
Finally, one interesting example where the modal behavior of a space was investigated also by means of wave-based acoustic simulation is the al Saflieni Hypogeum [11]. In this Neolithic structure, the authors investigated the possibility that the resonance frequencies observed in different chambers have been "tuned" in some way. Based on the on-site measurements, a numerical model was developed to test whether different dimensions of the chambers might have altered this tuning. However, as in the present study, this reconstruction reflects the "current" conditions which might have been different from the original conditions which, as the authors clearly state, may be difficult to imagine.

Conclusions
In this paper, the case of the rupestrian church of Saints Andrew and Procopius in Monopoli (Apulia, Italy) was considered. The church has a small volume and is in a very bad conservation state, but it was acoustically analyzed in its current state by means of measurements carried out according to ISO 3382-1 standard. A combination of different analysis methods, including not only simple room acoustical parameters, as in most previous studies, but also spectrograms, modal reverberation time, and detailed maps based on FDTD modelling, were used to interpret the results. Measurements pointed out a singular behavior characterized by short reverberation times (around 0.5 s) at frequencies above 250 Hz, and much longer values in the lowest bands reaching up to 2 s, also in combination with strong resonances due to the modal response of the room. The analysis of the narrow band spectra and of modal reverberation time confirmed such dependance and resonant frequencies, especially those appearing at 65 Hz and 79 Hz, which were explained as a function of the room dimensions (although with some approximations due to the strongly irregular shape). In order to acoustically simulate the space and analyze the distribution of sound both in the current and in any possible reconstructions of the original state, given the small dimensions that make it impossible to use geometrical acoustic models in the low-frequency range, the numerical FDTD method was applied by using a proprietary code developed in MATLAB. Thanks to a laser scanner survey, the geometry was simplified and voxelized so that FDTD could be applied. Results showed a very good agreement both in terms of predicted reverberation time and modal response, implying that FDTD method, despite its simplified implementation, is the most reliable approach to simulate such small spaces and obtain time-dependent responses. Consequently, a simulation of a possible reconstructed condition, characterized by plastered surfaces, debris removed from the floor, and some occupants in the space, was carried out. Results showed that the strong modal response still appeared, suggesting that liturgical singing might have taken advantage of such specific acoustic features.

Conflicts of Interest:
The authors declare no conflict of interest.