The paper introduces an innovative integration of multi-objective PSO and Pareto optimality for the joint inversion of MT and Rayleigh Wave Dispersion (RWD) data. The approach is rigorous, the synthetic tests are thorough, and the field application is well-motivated by the geotectonic context. However, several methodological, structural, and interpretative gaps (lacunae) diminish the clarity, reproducibility, and scientific robustness of the work. Addressing these issues will greatly enhance its publishability and reader understanding. Below are the comments to help the authors improve the work.

1. The paper references previous work using GA, cross-gradient, and petrophysical coupling. However, it does not quantitatively compare the proposed method with existing Pareto-GA techniques (e.g., Moorkamp et al. 2010), standard cross-gradient structural coupling, or traditional nonlinear inversion methods (Occam, Gauss–Newton). Without such benchmarks, it is difficult to evaluate the extent of the improvement the method offers. The authors could conduct at least one comparative test using a conventional inversion scheme to demonstrate computational or accuracy benefits.
2. In the synthetic results, especially Figures 6 and 8, seismic velocity is shown to be highly sensitive to noise, causing a shift of the Pareto front toward the MT axis. However, the work does not quantify noise thresholds at which model instability occurs, provide uncertainty analysis beyond posterior PDFs, or discuss robustness relative to the SNR typical in real MT/RWD surveys. It would be helpful if a noise-sensitivity analysis with multiple SNR levels could examine how stability depends on noise distribution.
3. The MT data reveal tensor skewness and non-1D features (page 17, Figure 10), and the geology is complex and heterogeneous. However, the paper relies solely on 1-D inversion for both MT and RWD. This could lead to misinterpretation because MT responses sometimes reflect 2D or 3D structures, and RWD curves across multiple paths naturally incorporate 2D structural heterogeneity. The joint inversion presumes a co-located subsurface structure, which is unlikely in tectonically complex regions. Please provide a more substantial justification for using 1-D modeling or include a preliminary 2D MT inversion or phase tensor analysis to demonstrate lateral homogeneity along the path.
4. The model uses 16 layers and 31 free parameters, but the MT and RWD datasets contain only about 20 periods each (page 11). This results in an underdetermined inversion. Although the constraint term Φc helps, there is no quantitative justification for selecting 16 layers, such as resolution analysis (e.g., sensitivity kernels, model covariance), or discussion of how many layers the data can practically resolve. Include a model resolution study or adopt an adaptive layering strategy, similar to an "Occam-style” smoothness approach.
5. The paper identifies Region A: Low resistivity and stable velocity as hydrothermal alteration; Region B: Low resistivity and low velocity as melt or fluids; and Region C: Low resistivity only as mineralization in the lower crust. These interpretations are plausible but not directly supported because no Vp/Vs, attenuation, or anisotropy data are used; no comparison with seismic tomography or MT 2-D sections is made; and no petrophysical temperature or pressure modeling is presented. Therefore, the interpretation might be overstated. If possible, include supporting geophysical evidence such as receiver functions, heat flow, or seismic tomography, or frame these interpretations more cautiously.
6. The study reports that one station requires three hours for 1,000 iterations (page 22), but it does not compare runtime with other algorithms, analyze scalability to 2D or 3D, or provide guidance on the number of particles versus computation time. These shortcomings indicate poor computational performance. Add a section on computational performance that discusses trade-offs, scaling predictions, and parallelization strategies.
7. There are no comparisons to borehole resistivity logs, local seismic velocity profiles, or previous MT/seismic surveys in the Biga Peninsula (many exist). This weakens confidence in the geophysical interpretation. Integrate or at least reference available independent constraints.
8. Some methodological steps lack clarity, such as how the initial particle positions were selected. How was convergence evaluated besides by iteration count? Why were c1 = c2 = 2.05 specifically chosen? How is the utopia point mathematically determined? These issues need to be addressed in the work.
9. Several figures (e.g., Pareto fronts, PDFs) lack explanations for colorbars, layer boundaries, and clear distinctions between the POS and mP* models. Error bars for observed data are missing, and there are inconsistencies in notation, such as ΦMT and ΦRWD appearing in different fonts. The terms “mP*,” “mP*-model,” and “PO-model” should be standardized. Layer numbering is inconsistent between the text and figures. Enhance figures to publication quality by improving annotations and standardizing color schemes.
10. Provide pseudocode for the algorithm.

This paper deals with the joint inversion of Magnetotelluric (MT) and Rayleigh wave dispersion (RWD) data for one‐dimensional resistivity- and velocity-depth models. In a new approach, the multiobjective particle swarm optimization (PSO) and the Pareto optimality method are utilized together, whereas in the past, global optimization algorithms with Pareto optimality were mostly integrated to jointly invert various geophysical data. The paper is theoretically clear in terms of a comprehensive description of the two applied methods, the creation of the synthetic dataset, and the references to previous applied joint inversion approaches. Despite some minor issues that need improvement, such as figure labels, inconsistent wording in the text, and mathematical symbols in the equations, I found the article quite straightforward to understand up to the description of the field datasets. One can already tell that the authors have established a valuable scientific basis for their application of joint inversion. However, the following chapters on one-dimensional joint inversion of MT and RWD inversion for a synthetic and real data example should be significantly improved in order to be accepted for publication. I suggest here a major revision of the paper and will present below my comments and questions that arose during my review and need to be revised.

Specific comments and questions:

MT:

Line 174-175: Why is the phase not considered as an input parameter for the inversion? MT is fundamentally based on the complex impedance tensor. The apparent resistivity is derived from the magnitude of the impedance tensor and the phase (derived from the argument of the complex impedance tensor) constrains the geometry and vertical resistivity gradients, which the apparent resistance alone cannot imply. By discarding the phase, valuable information about the depth localization of certain structures, such as conductors, will be lost and can lead to a non-uniqueness of the inverted resistivity-depth model. If this is the case, particularly due to heterogeneities near the surface, which usually lead to a static shift and are taken into account by the phase, the second dataset PWD must provide a strong structural constraint as compensation.

Line 174: The apparent resistivity is not a local or absolute value since it reflects a frequency-dependent, depth-averaged resistivity of the subsurface over the electromagnetic diffusion depth.

Line 209-210: Some things about the Electromagnetic (EM) fields have to be clarified here. Electric and magnetic fields are generated simultaneously by external EM sources, e.g. ionospheric/magnetospheric currents, as written in the paper. The propagating time-varying electric and magnetic fields interact with the conductive Earth, where they induce secondary currents within the Earth. The secondary currents can then be measured at the Earth’s surface.

Line 211: The MT period range has to be changed. The typical MT period range is 10^(-4) s (or 10 kHz) to 10^4 s.

Line 212: My comment in lines 209-210 indicates that EM waves are not only generated by magnetic fields. The sources of EM fields are electric charges and electric currents in the ionosphere or magnetosphere, such as lighting, substorms or solar winds. Sources are external to the Earth and internal currents are induced responses.

Line 213-216: The internal sources are not worth mentioning since they are not responsible for EM fields. The magnetization of a rock has nothing to do with MT diffusion processes.

Line 222-223: This sentence needs to be changed. The impedance tensor is not measured directly since the impedance is deduced from field measurements. As mentioned in my other comments before, the fundamental principle of the MT method is the estimation of the frequency-dependent impedance tensor from simultaneous measurements of the electric and magnetic fields at the Earth’s surface. The Earth modifies the ratio and phase between the electric and magnetic fields.

Line 223-224: How is the apparent resistivity derived from the complex impedance tensor? Perhaps the formula for apparent resistivity should be mentioned here, or at least that apparent resistivity is frequency dependent.

Line 226: Use the correct MT term “phase” and not “phase information”. (This should be set as a standard in the entire paper.)

Line 229-231: Consider a rewriting of the sentence to: “To estimate the impedance tensor, the time series were segmented into overlapping windows and transformed to the frequency domain using the FFT. Power and cross-spectra were then computed and stacked, from which the impedance tensor elements were estimated.”

Line 487-489: The authors have to mention at least once that for 1D MT modeling and inversion, a scalar or rotationally invariant impedance derived from the complete impedance tensor is used, since the impedance tensor itself is not used for modeling.

Line 484-485: As described in the paper, the MT data exhibit a 3D character beyond 1 s (not 10 s, as stated). The phase tensor ellipses become elliptical and change direction. In addition, the minimum phase angle changes abruptly in a short frequency/period range. Furthermore, the GULC station is located very close to the sea, which leads to a splitting of the apparent resistance curves of the two modes TE and TM at the longer periods of 10 s and above, here presumably Zxy and Zyx. This sea effect is caused by the strong conductivity contrast between seawater and the resistive crust/upper mantle, which can lead to a misinterpretation of the resistivities in the depth range of the lithosphere-asthenosphere boundary due to unreliable strong conductors in the resistivity-depth model. These periods must then be excluded for a 1D inversion. Therefore, the following questions must be answered more precisely in the paper:

1. How could the quality of the MT field data be improved by the ProcMT processing software? Are there many noise influences/effects visible in the recorded field data?
2. How do the main components of the impedance tensor look like? Are they close to zero, or large (hint to 3D MT data)?
3. Were the data rotated into a specific geoelectric strike direction after the processing to minimize the main diagonal elements of the impedance tensor and assume a 2D subsurface structure? Or were the MT field data treated directly as 1D (which would not be the accurate, at least for periods greater than 10 s, as can be seen from the phase tensor plots)?
4. Which component of the impedance tensor was ultimately used for the 1D inversion? Zxy or Zyx, since it is not mentioned?

General:

Line 255 & 282: The logarithm should be applied on the variable and not on the unit, such as log10(rho) [Ω𝑚] and log10[1, 5] Ω𝑚

Line 265: Why did the authors choose 31 model parameters with 16 layers for the models, even though the MT and RWD data sets only cover 20 periods? Normally, the layers should be set to a minimum, which also reduces the computing time, and then continuously increased after each iteration through adaptive layering. How many layers can the modeling algorithm resolve, which has not yet been shown?

Line 293: How are these values for c1, c2, k and chi defined? Are they derived from previous tests or inversions with the MOPSO and Pareto optimality method?

Line 319: How do the PDFs look for the others layers, e.g. in the resistivity-depth model? Since the PO-model shows greater deviations from the synthetic model, here at a depth of 7 or 10 km for the resistivity-depth model. And if the authors consider the 8th layer of the velocity-depth model to be incompatible, then the difference for the lowest resistivity at a depth of 10 km in the noise-free resistivity-depth model is already incompatible.

1. Have you tried the 1D inversion algorithm on a very simple three-layer model such as in Moorkamp et al. (2010) and verified the robustness and reliability of the algorithm? It appears that the inversion code has difficulty fitting the lowest resistivities, even in the case of field data that can be imaged most sensitively using the MT method. As mentioned earlier, this can happen when the phase of the MT data is not used in the inversion.
2. Can you additionally compare your joint inversion results with other 1D non-linear inversion codes for MT and RWD data? It would be beneficial to include an example in the paper that illustrates the accuracy of the 1D joint inversion models compared to other inversion methods using the same synthetic data.

Line 421-460: Here, the various geological structures of the Earth are identified using resistivity- and velocity-depth models. Region A can be attributed to volcanic rocks undergoing hydrothermal transformation. These findings are supported by earlier observations of the southeastern part of the Biga Peninsula, which were determined using seismic noise tomography. The other regions B, C, and < 5 km depth are explained by references to scientific studies of similar magmatic environments that could be responsible for mantle-derived melt, mineralization in the lower crust, and the sedimentary layer from the surface to a depth of 2 km on the Biga Peninsula. However, no connections are made to existing borehole logs or geophysical studies, such as Vp/Vs or S-wave tomography, previous MT surveys, and heat flow or fluid modeling studies that could support the statements about these different regions of the Earth. Some supporting studies of the research area would help to provide a more solid basis for the interpretations. In addition, the inclusion of the MT stations GURE and KULC in the geological map of Figure 4 could improve understanding of why more geothermal fields are found at the GURE site (in connection with the more permeable continental sediments). If possible, the location of the geothermal fields could also be mapped here.

Line 472-480: Instead of focusing the quantitative analysis of the sensitivity test on a single value, a representation of the actual data fit between the synthetic and modeled data (apparent resistivity and phase velocity) as shown in Figures 11 a) and b) would more clearly highlight the obvious changes to a better or worse data fit. A number for the data fit often says nothing about the actual influence of a substitution in the model parameters. This raises two further questions:

According to which criteria are the model parameters in the table selected?

What is the depth interval of the substituted model parameters, e.g., between 5 and 10 km depth for region A? (A visualization of the changes in the resistivity- and velocity-depth models would improve understanding about which sections of the depth models were changed for the sensitivity tests.)

Line 502-507: All these descriptions about the computing time of the PSO, the utilized computation infrastructure and the software must be mentioned in advance, for instance, directly in the “MOPSO and Pareto optimality parameters” section.

In general, a final chapter entitled “Conclusions” should be added to the paper, explaining and summarizing the main results and advantages of Pareto MOPSO in the joint modeling of MT and RWD data. In this context, lines 101–105 read like a concluding statement and can be included in the “Conclusions” section.

There is no section on “Data Availability” in the paper, or can the data only be published upon request?

Technical corrections:

Line 29: The abbreviations should be written out once in the main text: Magnetotelluric (MT) and Rayleigh wave dispersion (RWD)

Inconsistent wording (one has to be chosen for the entire paper): e.g. modelling to modeling; dataset or data set

Line 103: incorrect word

Line 107: Add PSO in the section title

Line 211: Maintain a consistent format for all units throughout the entire paper, here seconds or s.

Line 310: Should be 15%

Equation 1 & 6: The tensor product (Kronecker product) is not required here. Instead, element-wise multiplication with the Hadamard product must be inserted into both equations. Otherwise, the output would be a tensor that would increase the dimension, and the result would not be a vector in the model space.

Equation 3: The symbol for apparent resistivity in the equation does not match the one used in the following text.

Mathematical symbols are used in different notations in the paper, such as φ in line 140.

Figure 1: The search space is defined by an x-component on the x-axis. But what is on the y-axis? Shouldn't it be x1 on the x-axis and x2 on the y-axis?

Figure 2: The label b) is missing in the caption.

Figure 3: The apparent resistivity label in subfigure 3) looks distorted.

Figure 4: The color bar for the topographic map has no labels, e.g. “Elevation (m)”.

Figure 9: The labels a) and b) are either incorrectly assigned in the caption or in the figure itself.

Figure 10: The phase unit in the label should be specified in degrees and not as the phase symbol. Please use “phase” instead of “phase angle” here.

All subfigures are labeled (a), (b), etc. in the figure captions. However, the labels for the subfigures are a), b) and are also referred by these names in the text, which is inconsistent.