## 2 Computational Model  
A two-dimensional computational domain based on the experimental setup of [PERSON] and [PERSON] (2013); [PERSON] and [PERSON] (2013) is shown by the red-dashed box in Figure 2a. It covers a rectangular area of 0.132 m \(\times\) 7.4 m of the quarter-wavelength resonator. The mesh is designed to be denser within the vicinity of the heat exchanger (300 mm to the left and right from the heat exchangers) to resolve vortex structures issuing from and returning into the heat exchanger assembly. Elsewhere, a much coarser mesh is adopted as the problem consists of solving the governing equations for an empty resonator. The coarse mesh is illustrated in Figure 2b and the fine mesh in Figure 2c. This approach was validated against the theoretical predictions of linear thermoacoustic theory [PERSON] and [PERSON] (2013): 300 mm away from the heat exchanger, the pressure and velocity simply follow the linear model. A structured mesh of quadrilateral type was used for this model; the long computational domain is divided into several parts so that the distribution of mesh sizes could be made with acceptable quality. The minimum orthogonal quality was recorded to be 0.74 with maximum skewness of 0.31 (i.e., close to 1 for minimum orthogonality and close to 0 for maximum skewness [15]). The minimum mesh spacing in the \(x\)- and \(y\)-directions is 0.08 mm and 0.25 mm, respectively. The plates of the hot and cold heat exchangers (subsequently denoted as HHX for \"hot heat exchanger\" and CHX for \"cold heat exchanger\" throughout the paper) have thickness, \(d\), of 3 mm and are arranged in a parallel configuration with the spacing, \(D\), of 6 mm between them. This is shown in Figure 2d as a magnified view of a single channel between a pair of hot and cold plates; point'm' is the location of the joint where the hot and cold plates meet. The HHX and CHX are located next to each other with a joint positioned at 0.17\(\lambda\) from the pressure antinode. The wavelength, \(\lambda\), is defined as \(\lambda=a/f\), where \(a\) is the speed of sound and \(f\) is the frequency of the flow (13.1 Hz). Location of the pressure antinode for this 1/4-wavelength rig is at the end wall--cf. Figure 2a. The arrangement of the heat exchanger plates (10 pairs) follows the physical setup with five pairs of heated/cooled fins and five pairs of unheated ones, treated as adiabatic plates. The unheated plates provide consistent porosities for the flow (cf., [17]). The model extends over the whole height of the resonator to resolve the flow asymmetries observed experimentally [17]. The use of symmetry or periodic boundary conditions is inappropriate due to the presence of temperature-driven buoyancy effects. In addition, the experimental data are based on phases determined by comparing the phase of the velocity between the plates to the phase of the pressure at the pressure antinode. As the intention was to replicate as closely as possible the phenomena taking place in the physical apparatus, a full two-dimensional model was chosen. The boundary condition at the right end of the resonator is defined as a moving wall to replicate the acoustic displacement induced by the loudspeaker. The displacement, \(\delta\), is simply modelled as:  
\[\delta=\left[\frac{p_{a}}{\omega\rho_{m}a}\sin(k_{w}x_{s})\right]\cos(\omega t). \tag{1}\]  
Here, \(\omega\), \(\rho_{m}\), \(a\), \(k_{w}\), and \(t\) are the angular velocity, mean density, speed of sound, wave number, and time, respectively. The term \(p_{a}\) is the oscillating pressure at the location of the pressure antinode and \(x_{s}\) is the distance from the joint to the pressure antinode (shown as 4.6 m in Figure 2a). The mean pressure, \(p_{m}\), is 0.1 MPa. The drive ratio, \(p_{a}/p_{m}\), of the low amplitude flow investigated throughout this paper is 0.3%. The resonator wall was treated as adiabatic. The model was solved for two thermal conditions.