Dear Drs. Sarakorn &  Mukwachi,

I think this is a good, useful (and novel) method that you present. Also, the paper's well-written, well-organized, and illustrated well. I only have a small number of minor comments and suggestions.

Very last line of page 1, when talking about staggered-grid finite-difference: It would be good to be more precise in this statement. If it's the typical staggered-grid for the E-field equation, then it's the tangential electric field, I believe, that's continuous between cells, with no built-in condition on the current.

This is probably getting far too pedantic, so feel free to ignore. But, it would be possible to formulate the typical integral equation approach (whether one that doesn't work for high contrast or one that does) using, e.g., unstructured tetrahedral meshes, to accurately model a real-life ore-body, etc. And it is possible to formulate staggered-grid finite-difference (or "finite-volume") for, e.g., unstructured tetrahedral meshes. I'd argue, therefore, that it's not so much the methods themselves that are limited to rectangular discretizations but the implementations that people have derived and coded up.

Lines 37-38: One problem with the simple rectilinear meshes is the poor mesh quality, and hence ill-conditioning of the system, resulting from extending padding zones out to the boundary of the domain, which ends up with very pizza-box-like cells. This is very nicely illustrated by your Figure 5.

What field do you use for your boundary conditions: the electric field in a homogeneous halfspace (and if so, what's the conductivity of the halfspace?), or in a horizontally layered Earth, or in 2D models?

Why use a mesh made from quadrilateral elements and not one made from triangular elements, and then extrude those triangles downwards?

Lines 96-97: This is perhaps slightly important. Can your meshing, and hence the forward-modelling code, handle cells of variable vertical extent? This is what the sentence over these two lines suggests. However, there is nothing else in the manuscript that mentions this. If your method (meshing & MT modelling) can handle this, it would be good to see this shown off and demonstrated in an example (perhaps the trapezoidal hill model of Nam et al.).

Figure 1: In this figure you're meaning to indicate on which boundaries of the domain the (non-zero) tangential component of the background electric field is being applied for the two different polarizations, is that correct? But you still need to provide conditions on the other boundaries, it's just that these involve forcing the tangential electric field to zero?

Line 149: Can't you use analytic formulae for the integrals? Why use a numerical integration technique? Doesn't this take quite a bit longer?

Line 165 and thereabouts: It's not 100% clear: Do you use a sparse, direct solver from NumPy (or elsewhere), or an iterative solver? (If an iterative solver, are you doing any kind of divergence correction, or using a specially designed preconditioner, to help with convergence?)

Figure 5: What does the vertical cross section at y=0km look like? That would be interesting to see, including the cells up under the observation locations. Presumably the aspect ratios of those cells are good? Same for Figure 7.

Figure 11: It would be interesting to see a vertical slice through the mesh along the line of observation locations. Do the vertical extents of the cells in a layer vary?

Figure 12: I think it would be good to re-do this figure. Perhaps digitize the data points from the figure in Nam et al. (their Figure 12, isn't it?) as this would allow you to plot their data on a new graph that you then have total control over. And don't include the 2D results, which are the ones that have the elevated apparent resistivities in the middle of the survey line for the yx-polarization, counter to what the 3D results are doing: this is just distracting. Even having the crosses (2D results) on the graph of xy apparent resistivities is distracting and takes away from how close your results are to the 3D results of Nam et al.

Best wishes,

Colin Farquharson.

Dear Drs. Sarakorn and Mukwachi,

This is a well-written and thorough manuscript that makes a valuable contribution to the field of 3D MT modelling. I only have a few of minor comments and suggestions.

From the presented material it is difficult to intuitively understand how the semi-unstructured horizonal mesh relates to the more conventional vertical mesh, especially in the first few model layers. I suggest that it would be very instructive for a reader to have an expanded version of Figure 11. The expanded figure would include some additional horizontal slices through the shallow parts of the model showing the mesh structure. Alternatively, a slice through the model under the topography feature with a depth extent of approx. 4 km would allow a fuller understanding of how the mesh functions in the shallow subsurface.

A similar additional section or couple of slices might also be instructive for the mesh structure in Figure 7.

Additional minor suggestions -

There are a few overly long sentences in the manuscript that could be shortened for clarity. Specifically:

• Line 21. “In the early era, the IE method, based on applying Green’s function to the scattering equation, and the SGFD, a variant of the SGFD method that enforces the continuity of electric current along the block’s edges and across its faces (Yee, 1966; Siripunvaraporn et al., 2002), are accurate and efficient for the simple 3D domain.“
• Line 94. “By combining the geometric flexibility of a paving algorithm for the horizontal plane with controlled, adaptable vertical layering, this approach is expected to offer new and improved performance for 3D MT modeling through the EBFE approach, providing more accurate solutions and better computational efficiency by effectively addressing electromagnetic fields across highly contrast geophysical structures, which are crucial for correct simulation of MT responses.”
• Line 104. “By combining the geometric flexibility of a paving algorithm for the horizontal plane with controlled, adaptable vertical layering, this approach is expected to offer new and improved performance for 3D MT modeling through the EBFE approach, providing more accurate solutions and better computational efficiency by effectively addressing electromagnetic fields across highly contrast geophysical structures, which are crucial for correct simulation of MT responses.”

Other comments:

• The sentence on line 160 needs to be edited to make the meaning clear.
• The sentence on line 229 needs to be edited to make the meaning clearer.
• Nam et al. is incorrectly marked as Name et al. in Figure 12.
• Figure 5 uses both east-west/north-south and x/y nomenclature when describing the model in the bottom two panels. It would be better to only use one reference system to make the figure clearer.

All the best.