Heavy-duty land-based gas turbines with can-annular combustors typically exhibit a set of thermoacoustic modes with similar frequencies, due to their rotational symmetry and the weak coupling between neighbouring cans. The Bloch boundary condition is often used to simplify these combustors. Because the widely used (quasi-)one-dimensional model cannot naturally consider high-order modal components in the three-dimensional acoustic field, it is valid only when the frequency is sufficiently low and the length of the can-to-can interaction gap is sufficiently short. In this paper, we present a theoretical model using the mode-matching method to consider high-order modal components. A reflection coefficient matrix is used to accurately describe the acoustic characteristics of the gap area. This reflection coefficient matrix is found to depend on the nondimensional cross-sectional area ratio, the Helmholtz number based on the length of the gap area, and the Bloch number. To facilitate a full low-order network modelling of the high-frequency thermoacoustic modes in can-annular combustors, a flame transfer matrix is also used to describe the modal coupling effect across the flame. Based on this model, the thermoacoustic problem of a four-flame can-annular combustor is studied as an example. It is shown that the frequency and growth rate of the thermoacoustic modes can be predicted accurately using the reflection coefficient matrix and flame transfer matrix. The effects of high-order modal components across the flame and in the gap areas are studied separately. The results show that the model containing both of the two coupling effects has the widest range of validity. These couplings can be ignored if only low-frequency plane acoustic waves are considered but need to be carefully considered for higher frequency modes. Moreover, when multiple flames are considered, asymmetric flame distributions are found to more easily lead to modal coupling compared to symmetric flame distributions. All results are validated by comparing with the numerical solutions of the full two-dimensional Helmholtz equation using the finite element method.