###### Abstract  
Most geophysical inversions face the problem of non-uniqueness, which poses a challenge in the mapping and delineation of the subsurface anomalies. To tackle this challenge, a combined local and global optimization approach is considered for jointly inverting two-dimensional direct current resistivity (DCR) and seismic refraction (SR) data that aim to estimate the corresponding physical model parameters. In this combined approach, the output of the local optimization method is used to determine the search space and tuning parameters for the global optimization algorithm. The multi-objective genetic algorithm (non-dominated sorting genetic algorithm) was utilized to jointly optimize the objective functions of two different methods. Because the genetic algorithm is a population-based optimization method, it requires numerous forward calculations. To deal with the expected high computational cost associated with this approach, parallel computing was utilized for the forward function evaluations to reduce the run time of the entire process. The proposed approach was tested using synthetic two-dimensional resistivity and velocity models that had three different types of anomalies (dyke, positive, and combined positive and negative). The results showed an improvement in the anomaly delineation in the output of the combined local and global optimization method compared with the local optimization method. Additionally, similar synthetic models were tested using only the single objective global optimization algorithm (conventional global optimization), which showed promising anomaly delineation. 1 (2022)  
function formed by a weighted or arithmetic sum of the individual objective functions of corresponding methods [23]. However, the issue of parametric coupling (i.e., model integration) arises when we are dealing with the inversion of datasets acquired with more than one geophysical method [24]. In this study, we apply the structural model coupling approach, which involves the cross-gradient constraint method that is commonly used for geophysical inversion [8]. The basic idea using this approach is that the gradients of relevant model parameters are spatially correlated. The cross-gradient approach has been adapted for the joint inversion of different geophysical data [13; 14; 15; 16; 25; 26; 27; 28; 29]. [PERSON] et al. [17] applied the cross-gradient algorithm in a joint inversion involving both the controlled-source audio magnetotelluric (CSAMT) and magnetic methods. [PERSON] et al. [15] formulated an objective function concerning weighted cross-gradient, which limits the dominance of one type of model parameter to another. [PERSON] et al. [30] utilized the cross-gradient constraint to impose a common structural framework in the joint inversion of EM and acoustic data, which reconstructs the structures satisfactorily. Furthermore, [PERSON] et al. [31] applied a cross-gradient technique to invert magnetotelluric (MT) and gravity data and tested their algorithm using both synthetic and real datasets. Finally, [PERSON] et al. [32] introduced a new approach to the cross-gradient constraint, which involved the use of an irregular grid in unstructured mesh in the finite element method. They used the method to invert the DCR and ground penetrating (GPR) data. All the above examples used the conventional or local optimization method that incorporated all acquired data within the same Jacobian using a unique objective function for all geophysical methods. The local optimization techniques are applied iteratively to obtain an updated model that minimizes the objective function, which may not be the global solution to the inverse problem. The global optimization algorithm in geophysical inversion might be used to search for a solution space to avoid being stuck in the local minimum of the objective function. For instance, [PERSON] et al. [33] applied the particle swarm optimization (PSO) algorithm in a parametric inversion involving magnetic data. [PERSON] et al.