Determination of Optimal Diffusion Coefﬁcients in Lake Zirahu é n through a Local Inverse Problem

: In this study, optimal diffusion coefﬁcients for Lake Zirahu é n, Mexico, were found under particular conditions based on images taken with a drone of a dye release experiment. First, the dye patch concentration was discretized using image processing tools, and it was then approximated by an ellipse, ﬁnding the optimal major and minor axes. The inverse problem was implemented by comparing these observational data with the concentration obtained numerically from the 2D advection–diffusion equation, varying the diffusion tensor. When the tensor was isotropic, values of K 11 = K 22 ≈ 0.003 m 2 /s were found; when nonequal coefﬁcients were considered, it was found that K 11 ≈ 0.005 m 2 /s and K 22 ≈ 0.002 m 2 /s, and the cross-term K 12 inﬂuenced the results of the orientation of the ellipse. It is important to mention that, with this simple technique, the parameter estimation had consequences of great importance as the value for the diffusion coefﬁcient was bounded signiﬁcantly under particular conditions for this site of study.


Introduction
Lakes are vitally important components that provide essential ecosystem services, such as water for drinking, and food supply and sites for recreation and tourism [1]. However, the quality of lakes' water resources is deteriorating rapidly due to the load or spread of pollutants, excessive nutrient inputs, and all types of sediments. In particular, suspended solids affect dissolved oxygen levels and the temperature, interfering with the mixing and decreasing the dispersion processes toward deeper layers [2][3][4]. The sediment transport estimation is often difficult, time-consuming, and expensive [5]. Nonetheless, observations and monitoring are important aspects of ecosystem services that can be used to map the distribution of sediments and can also be used as a target for numerical models to predict sediment movement and spreading, in particular, knowing that the diffusion coefficient is essential for the numerical modeling of such transport processes [6][7][8][9].
The transport of pollutants can be described by many factors, including advection, diffusion, dispersion, reaction, dilution or mixing, retardation, and decay. These factors are usually incorporated into transport equations, which describe phenomena such as mass transfer or heat, fluid, waves, and momentum transfer [10,11]. Nevertheless, this equation is usually solved numerically for advection-diffusion terms solely, and such numerical models are frequently implemented with known constant coefficients [10,11].
The parameterizations of a set of coefficients have been the topic of ongoing research, most of it carried out empirically and, more recently, based on theoretical considerations with the use of the inversed problem approach. The use of adequate formulations in inverse modeling has proven to be a highly useful method for estimating these parameters and improving the fit of observational data [11][12][13][14][15]. With the use of parameter estimation for a specified set of data, the numerical modeling of the transport equation can be improved for a particular study site.
Likewise, dye tracers are a highly useful tool to estimate advection-diffusion characteristics implemented to measure several aspects of the transport process, such as to measure certain effluent discharge rates and horizontal dispersion coefficients; it is also used to study circulation patterns and to calibrate-validate numerical models used in the forecast activity [10,[16][17][18][19]. Along with dye-tracers, aircraft systems (drones) and image processing tools can be used to visualize the dye plume and provide accurate tools to identify, map, and predict the spread of hazardous agents [18,20,21].
For image processing or pattern recognition, several algorithms have been used to detect objects fitting certain geometrical shapes, such as circles or ellipses, to the images. Although there are several sophisticated algorithms for detecting shapes from images, the least squares method is the most often present in these algorithms [22,23].
In this work, a set of images of a dye tracer patch was used to determine the optimal diffusion coefficients from the advection-diffusion equation. The hypothesis was that, based on the inverse problem approach and the observational data of a dye patch, as a proxy for sediment transport, we expected to find optimal diffusion coefficients for the study site and under particular conditions without the involvement of sophisticated measurements to estimate it.

Materials and Methods
To generate data for solving the advection-diffusion equation, a dye tracer experiment was performed in Lake Zirahuén, Mexico. In general, the methodology applied in this study consists of (i) the estimation of the dye distribution using image processing tools to obtain the area of the dye patch based on an optimization technique, and (ii) the solution of the nonlinear inverse problem for the advection-diffusion equation in order to find the optimal diffusion coefficients by comparing this analytical solution against the observational dye patch distribution; see Figure 1 for an overview and a better understanding of the study. for a specified set of data, the numerical modeling of the transport equation can be improved for a particular study site. Likewise, dye tracers are a highly useful tool to estimate advection-diffusion characteristics implemented to measure several aspects of the transport process, such as to measure certain effluent discharge rates and horizontal dispersion coefficients; it is also used to study circulation patterns and to calibrate-validate numerical models used in the forecast activity [10,[16][17][18][19]. Along with dye-tracers, aircraft systems (drones) and image processing tools can be used to visualize the dye plume and provide accurate tools to identify, map, and predict the spread of hazardous agents [18,20,21].
For image processing or pattern recognition, several algorithms have been used to detect objects fitting certain geometrical shapes, such as circles or ellipses, to the images. Although there are several sophisticated algorithms for detecting shapes from images, the least squares method is the most often present in these algorithms [22,23].
In this work, a set of images of a dye tracer patch was used to determine the optimal diffusion coefficients from the advection-diffusion equation. The hypothesis was that, based on the inverse problem approach and the observational data of a dye patch, as a proxy for sediment transport, we expected to find optimal diffusion coefficients for the study site and under particular conditions without the involvement of sophisticated measurements to estimate it.

Materials and Methods
To generate data for solving the advection-diffusion equation, a dye tracer experiment was performed in Lake Zirahuén, Mexico. In general, the methodology applied in this study consists of (i) the estimation of the dye distribution using image processing tools to obtain the area of the dye patch based on an optimization technique, and (ii) the solution of the nonlinear inverse problem for the advection-diffusion equation in order to find the optimal diffusion coefficients by comparing this analytical solution against the observational dye patch distribution; see Figure 1 for an overview and a better understanding of the study.

The Inverse Problem
To solve the inverse problem, the nonlinear least squares method was used. These problems arise in the context of fitting a parameterized mathematical model to a dataset to minimize an objective function that can be expressed as follows: where r i is a smooth function from R n to R. Such a minimization problem is derived from curve fitting by least squares, where f (x) measures the difference between the model and the observational data, i.e., r i represents the residuals.
Although there are many optimization methods, most of them have the same structure based on Newton's method. To solve Equation (1), the Levenberg-Marquardt method that combines two algorithms, the gradient descent method and Gauss-Newton method, was implemented. The Levenberg-Marquardt method is more similar to a gradient descent method when the parameter is far from the optimal and works as a Gauss-Newton method when the parameters are close to the optimal [24-28].

Image Processing and the Inverse Problem to Estimate the Concentration
During a field campaign conducted in July 2018 in Lake Zirahuén, a dye release experiment and an unmanned aerial vehicle (DJI Phantom 4 drone) were used to capture videos and images of a dye plume. Images were captured from 50 m high at a rate of approximately one image every 3 s. In this experiment, natural food colorant (not harmful for the environment) was used as a tracer. The method consisted of the instantaneous release of the dye on the surface used to determine the transport.
Image processing was performed by adjusting the red, green, and blue (RGB) values to increase the contrast of every image and select only the dye plume. Then, working with the binary mask-image and using an edge detection approach, the boundary points of the region were determined to evaluate the extension of the patch. The dye patch area was estimated based on the hypothesis that, after the release, the dye cloud can be approximated by an ellipse, as the natural propagation of the tracer is better characterized by this geometrical form.
To estimate the proposed concentration ellipse (C e ) that fits the dye tracer patch, it was first necessary to find the center and orientation of the data considering the algebraic equation of the general conic to define the minimization problem. In this case, the nonlinear optimization problem was solved for finding the optimal lengths of the major and minor axes.
To obtain the ellipse parameters, we considered only the boundary points, (ζ i , η i ), i = 1, 2, . . . n, of the region. Then, if the geometrical center (ζ c , η c ) is defined as: Then, based on the variance: and the covariance: the orientation of the principal axes is determined by: Mathematics 2021, 9, 1695 4 of 14 and the variance in the major and minor axes by: To calculate the optimal length of the major axis (x 1 ·σ ma ) and the minor axis (x 2 ·σ mi ), the ellipse equation considered with the center in (ζ c , η c ) and θ, the rotation angle: was used to define the quadratic error of Equation (1) for all the points, where r i (x 1 , x 2 ) is the residual defined by (7) when (ζ, η) is replaced with (ζ i , η i ). Therefore, this equation was minimized based on the two parameters (x 1 and x 2 ) to obtain the best ellipse that represents the dye tracer area [23,29,30].

Mathematical Formulation of the Direct and Inverse Problem of the Transport Equation
The mathematical formulation most suitable for simulating the phenomena at hand of the transport of a substance in a fluid is the advection-diffusion equation: with initial and boundary conditions: where ϕ represents the tracer concentration and n is an outward vector normal to the boundary of the spatial domain Ω; v is the bidimensional velocity field; K is the diffusion tensor; X = (ζ, η) and t are the spatial and temporal variables, respectively. The forward problem for the advection-diffusion equation was solved using the MATLAB PDE Toolbox (The MathWorks, Inc., Natick, MA, USA). The core of this toolbox algorithm uses the finite element method for problems defined on bounded domains. Here, the advection term was provided as a function and was set as the f-coefficient of the PDE Toolbox.
The general procedure to solve this equation is to choose the flow field and a diffusion coefficient in advance and compute the tracer distribution by solving the forward problem; then, the nonlinear inverse problem is solved to obtain the optimal diffusion coefficients based on the observational data from the spatial distribution of the tracer and the flow field.
Knowing the diffusion coefficient that is essential for internal processes, such as stratification and mixing [6][7][8][9], can be achieved from the diagonal tensor for diffusion, expressed as: and if the flow is not parallel to the cartesian axes, then the cross-symmetric terms are required: To find the optimal diffusion coefficients for Lake Zirahuén, the inverse problem consists of searching the diffusion coefficients that yield the best possible fit to the observed values by solving the problem in the sense of a nonlinear least squares problem [11,[31][32][33], as in Equation (1), where: r = ϕ − ϕ obs ,K = K 11 = K 22 and K 12 = K 21 = 0or K 11 = K 22 and K 12 = K 21 = 0. (12) This equation makes an optimal approximation by comparing the value of concentration ϕ obtained by solving the advection-diffusion Equation (8), with the observational data ϕ obs . The model is rerun with the the updated K value until the minimal difference between the simulated (ϕ) and observed (ϕ obs ) is found.
For practical solutions of the minimization of the objective functional Equation (1), the MATLAB subroutine lsqnonlin is used, starting at an initial estimate based on the trust-region-reflective algorithm or the Levenberg-Marquardt method [34][35][36][37][38]. A good convergence indicator is given by a small first-order optimality measure and a moderate number of iterations.

The Velocity Field Calculation
The measurements to obtain the flow velocities are usually made by sensors, such as a current meter and acoustic profilers, but this equipment is expensive, which could be an inconvenience [29]. An alternative method is the use of unmanned air vehicles (UAVs) that allow high-resolution data collection in larger areas at minimal cost. This technology was used to measure the velocity field on the surface by using image-based techniques, resulting in a cost-effective and environmentally friendly method [21,39].
Then, from the contours of the dye patch, the mean velocity field was also calculated, to solve the 2-D advection-diffusion equation. The velocity field was obtained by comparing two consecutive patch contours, considering the geometrical center of the ellipses (x c , y c ) with a time step of 40 s. From the geometrical center of the ellipses, the distance between consequent contours was obtained, and then the velocity was computed by dividing the distance for the time step between contours. In this way, velocity components were calculated as follows: and the average velocity for each component would be:

The Setting of the Model
To solve the direct problem of the transport equation, from the information of the dye tracer experiment, the region Ω = {(ζ, η) : 0 ≤ ζ ≤ 30, 0 ≤ η ≤ 17} was considered according to the number of pixels of the image and the height from where the images were captured. The velocity field (u, v) and the diffusion coefficient K were considered constants. The Neumann boundary condition was set to zero.
The initial condition was the area related to the ellipse of the first image set as: where σ 2 ma = 0.255 m 2 ; σ 2 mi = 0.498 m 2 ; θ = 0.4876 rad; ζ c = 16.48 m; η c = 6.47 m. The domain was portioned into unstructured triangular elements with 2447 nodes and 1176 elements, for elements of approximately 1 m, as shown in Figure 2. For the solver options, the absolute, relative, and residual tolerance were set as the default values. The accuracy of the simulation was checked by reducing the relative tolerance, and with simulations with different element sizes (not shown).
accuracy of the simulation was checked by reducing the relative tolerance, and with simulations with different element sizes (not shown).

Image Processing
In Figure 3a, the scheme of the image treatment, which consists of the selection of (i) the RGB (red, green, and blue) intensity, (ii) the inverted mask image, (iii) the binary image, and (iv) the selection only of the dye patch, is shown. After georeferencing the image and calculating the geometrical center and the optimal axes, in Figure 3b, the binary image with the boundary dataset of the dye tracer patch, along with the optimal ellipse fit, is shown. In Figure 4, the main characteristics of the dye patch extent mapping are shown. The images show a dye patch dispersing over time. The first image shows the dye patch after several seconds of release, whereas the last was taken before the patch was not identifiable, due to the reflection of the cloudy sky. The evolution of the corresponding binary image and the optimal ellipse is shown on the right side of Figure 4a-h. It is worth noting that during this experiment (19 July 2018), a light breeze was registered of less than 4 km/h, but enough to displace the boat that was drifting freely with the motor turned off so as to not interfere with the dye patch. Other conditions measured during this experiment were the air temperature of approximately 20 °C, the surface water temperature of

Image Processing
In Figure 3a, the scheme of the image treatment, which consists of the selection of (i) the RGB (red, green, and blue) intensity, (ii) the inverted mask image, (iii) the binary image, and (iv) the selection only of the dye patch, is shown. After georeferencing the image and calculating the geometrical center and the optimal axes, in Figure 3b, the binary image with the boundary dataset of the dye tracer patch, along with the optimal ellipse fit, is shown.

Image Processing
In Figure 3a, the scheme of the image treatment, which consists of the selection of (i) the RGB (red, green, and blue) intensity, (ii) the inverted mask image, (iii) the binary image, and (iv) the selection only of the dye patch, is shown. After georeferencing the image and calculating the geometrical center and the optimal axes, in Figure 3b, the binary image with the boundary dataset of the dye tracer patch, along with the optimal ellipse fit, is shown. In Figure 4, the main characteristics of the dye patch extent mapping are shown. The images show a dye patch dispersing over time. The first image shows the dye patch after several seconds of release, whereas the last was taken before the patch was not identifiable, due to the reflection of the cloudy sky. The evolution of the corresponding binary image and the optimal ellipse is shown on the right side of Figure 4a-h. It is worth noting that during this experiment (19 July 2018), a light breeze was registered of less than 4 km/h, but enough to displace the boat that was drifting freely with the motor turned off so as to not interfere with the dye patch. Other conditions measured during this experiment were the air temperature of approximately 20 °C, the surface water temperature of In Figure 4, the main characteristics of the dye patch extent mapping are shown. The images show a dye patch dispersing over time. The first image shows the dye patch after several seconds of release, whereas the last was taken before the patch was not identifiable, due to the reflection of the cloudy sky. The evolution of the corresponding binary image and the optimal ellipse is shown on the right side of Figure 4a-h. It is worth noting that during this experiment (19 July 2018), a light breeze was registered of less than 4 km/h, but enough to displace the boat that was drifting freely with the motor turned off so as to not interfere with the dye patch. Other conditions measured during this experiment were the air temperature of approximately 20 • C, the surface water temperature of 21 • C, and a stratified water column, creating almost ideal conditions to carry out the experiment [40]. 21 °C, and a stratified water column, creating almost ideal conditions to carry out the experiment [40].

The Velocity Field
The spatial evolution of the dye patch shows a tendency to the right due to the wind forcing at the site. It was observed that the dye patch was advected ~6 m with an almost homogeneous dispersion of approximately 1 to 5 m of extension for 290 s. The positions of the geometrical center at each time are shown in Table 1. From the comparison of two consecutive patch contours, the velocity data were obtained according to Equation (14); see Figure 5.

The Velocity Field
The spatial evolution of the dye patch shows a tendency to the right due to the wind forcing at the site. It was observed that the dye patch was advected~6 m with an almost homogeneous dispersion of approximately 1 to 5 m of extension for 290 s. The positions of the geometrical center at each time are shown in Table 1. From the comparison of two consecutive patch contours, the velocity data were obtained according to Equation (14); see Figure 5.

General Aspects
The primary goal of the inverse problem modeling approach is to determine the parameter values that yield the best fit to the observed data. In Figure 6, the objective function ( ) of Equation (1) and the number of iterations are presented. A smooth convergence is shown, in which the minimization process is completed with a prescribed tolerance of = 1 × 10 −6 . From the analyzed cases, it was indicated that, when 11 ≠ 22 (Table 2), the computational time was greater and the number of iterations was 1 or 3 times higher, but the solution was more precise as the error decreased. When 11 = 22 , the solution was obtained after seven iterations. Additionally, when the anisotropic terms were considered, when 11 ≠ 22 and 12 = 21 = 0, the convergence was achieved after eight iterations, and when 11 ≠ 22 and 12 = 21 ≠ 0, the optimal solution was obtained after ten iterations.

General Aspects
The primary goal of the inverse problem modeling approach is to determine the parameter values that yield the best fit to the observed data. In Figure 6, the objective function f (K) of Equation (1) and the number of iterations are presented. A smooth convergence is shown, in which the minimization process is completed with a prescribed tolerance of tol = 1 × 10 −6 .

General Aspects
The primary goal of the inverse problem modeling approach is to determine the parameter values that yield the best fit to the observed data. In Figure 6, the objective function ( ) of Equation (1) and the number of iterations are presented. A smooth convergence is shown, in which the minimization process is completed with a prescribed tolerance of = 1 × 10 −6 . From the analyzed cases, it was indicated that, when 11 ≠ 22 (Table 2), the computational time was greater and the number of iterations was 1 or 3 times higher, but the solution was more precise as the error decreased. When 11 = 22 , the solution was obtained after seven iterations. Additionally, when the anisotropic terms were considered, when 11 ≠ 22 and 12 = 21 = 0, the convergence was achieved after eight iterations, and when 11 ≠ 22 and 12 = 21 ≠ 0, the optimal solution was obtained after ten iterations.  From the analyzed cases, it was indicated that, when K 11 = K 22 (Table 2), the computational time was greater and the number of iterations was 1 or 3 times higher, but the solution was more precise as the error decreased. When K 11 = K 22 , the solution was obtained after seven iterations. Additionally, when the anisotropic terms were considered, when K 11 = K 22 and K 12 = K 21 = 0, the convergence was achieved after eight iterations, and when K 11 = K 22 and K 12 = K 21 = 0, the optimal solution was obtained after ten iterations.
In Figure 7, the time evolution of the case when K 11 = K 22 and K 12 = K 21 = 0 is shown as an example; the image corresponds to the direct problem for the optimal diffusion coefficients found. It is shown that the evolution of the concentration moved to the right due to the given velocity field, and as time progressed, a lower concentration was observed in the center.   Figure 7, the time evolution of the case when 11 ≠ 22 and 12 = 21 = 0 is shown as an example; the image corresponds to the direct problem for the optimal diffusion coefficients found. It is shown that the evolution of the concentration moved to the right due to the given velocity field, and as time progressed, a lower concentration was observed in the center.

The Final Concentration Snapshots
The distribution of the concentration for the final step solution from the simulation is shown in Figure 8a in color contours. In order to obtain a simpler visual analysis, one arbitrary contour is displayed by a black line for the observations and by a red line for the numerical solution of the transport equation once the optimal coefficient was found.
The distribution of the concentration when 11 = 22 had a circular shape (as expected) and moved slightly to the right with respect to the center of the observations (the ellipse) (Figure 8a). The optimal diffusion coefficients in this case were 11 = 22 ≈ 0.003 m 2 /s. Figure 7. Concentration distribution. Concentration contours obtained with the optimal diffusion coefficient K 11 = K 22 and K 12 = K 21 = 0. The gray contour exemplifies the temporal evolution.

The Final Concentration Snapshots
The distribution of the concentration for the final step solution from the simulation is shown in Figure 8a in color contours. In order to obtain a simpler visual analysis, one arbitrary contour is displayed by a black line for the observations and by a red line for the numerical solution of the transport equation once the optimal coefficient was found. When 11 ≠ 22 and 12 = 21 = 0, the solution improved as the resulting numerical shape was an ellipse similar to the one proposed in the observations (Figure 9a), but slightly shifted to the right. The distribution of the concentration when K 11 = K 22 had a circular shape (as expected) and moved slightly to the right with respect to the center of the observations (the ellipse) (Figure 8a). The optimal diffusion coefficients in this case were K 11 = K 22 ≈ 0.003 m 2 /s. When K 11 = K 22 and K 12 = K 21 = 0, the solution improved as the resulting numerical shape was an ellipse similar to the one proposed in the observations (Figure 9a), but slightly shifted to the right. When 11 ≠ 22 and 12 = 21 = 0, the solution improved as the resulting numerical shape was an ellipse similar to the one proposed in the observations (Figure 9a), but slightly shifted to the right. Additionally, in Figure 10a, when 11 ≠ 22 and 12 = 21 ≠ 0, a major improvement was achieved as the solution had an ellipse shape, and the orientation matched with the observations; however, a slight shift to the right was still presented. Additionally, in Figure 10a, when K 11 = K 22 and K 12 = K 21 = 0, a major improvement was achieved as the solution had an ellipse shape, and the orientation matched with the observations; however, a slight shift to the right was still presented. In Figures 8b, 9b, and 10b, the error between the data and the fit, characterized by the sum of square errors (SSE) and the determination coefficient ( 2 ), is shown. In Figure 8b, SSE = 0.93 and 2 = 084; meanwhile, in Figure 9b, the values correspond to SSE = 0.71 and In Figure 8b, Figure 9b, and Figure 10b, the error between the data and the fit, characterized by the sum of square errors (SSE) and the determination coefficient (R 2 ), is shown. In Figure 8b, SSE = 0.93 and R 2 = 084; meanwhile, in Figure 9b, the values correspond to SSE = 0.71 and R 2 = 0.88, and in Figure 10b, SSE = 0.67 and R 2 = 0.89. Then, the best result was quantitatively obtained for the case when K 11 = K 22 and K 12 = K 21 = 0, as R 2 was closer to unity, and the error tended to zero ( R 2 → 1 y SSE → 0 ).

Discussion
Based on the results of the diffusion coefficients, the value for these parameters in the study area can be significantly bounded, which is crucial for the forward and inverse problems, as in several studies, the considered values are in the range between 0.0001 and 1 m 2 /s, and such a wide range can be reflected in significant changes in the dynamics of the transport phenomena. For example, studies [41][42][43][44][45] have been carried out where the distribution of pollutants in lakes and ponds is analyzed using numerical models that solve the advection-diffusion equation, and different values for the diffusion coefficient are considered that span several orders of magnitude (see Table 3). Table 3. Studies with different values for the diffusion coefficients considered.

References Diffusion Coefficient
Kusuma et al. It should be noted that in some of these studies, information was obtained from data of the velocity field, as well as the concentration of different substances measured with very sophisticated instruments, and in several studies, they considered that the propagation rate is proportional to the concentration gradient to estimate the diffusion coefficient.
In the present study, it should be noted that a more dynamic approach was used in which the transport equation itself was basically used to determine its corresponding and more adequate value for a set of observations based on an inverse problem approach.

Conclusions
From the dye tracer images, it was noticed that the shape of the distribution of the tracer was not uniform, but for numerical purposes, the ellipse resulted in an improved and accurate approximation; as the typical consideration was a circular shape, although showing similar results, the ellipse was better. Additionally, a perfect match with the ellipse approximation was not expected, as there would be more physical factors to consider in reality, for example, if there was no wind, the diffusion of a point source in a uniform medium should be perfectly circular; however, there are other physical mechanisms that make the diffusion parameter variable according to the conditions such as turbulence, enhanced by waves and its secondary circulations (e.g., Stokes drift and Langmuir circulations). Other factors that could be considered are the chemical or biological properties in the waterbody; however, with further development and measurements, this work can be improved.
Nonetheless, the result of the simulations when K 11 = K 22 was K ≈ 0.003 m 2 /s, whereas, when nonequal coefficients were considered, K 11 = K 22 , K 11 ≈ 0.005 m 2 /s, and K 22 ≈ 0.002 m 2 /s were found. The value of the coefficients showed no significant changes, regardless of whether the term K 12 (= K 21 ) was considered; however, this parameter was important as it influences the orientation of the ellipse.
Counting with measurements is of great importance as the results of different studies have indicated that pollutants move with the flow and diffuse over time. However, sometimes, the resources are not available to carry out these measurements, hence the importance of this study in terms of the sampling technique and the parameter estimation, which were relatively simple to use. In addition, the values obtained in this paper contribute to a better understanding of the transport phenomena in Lake Zirahuén, which have not been extensively analyzed to date.
Even though the aim was to have a practical implementation to find and bound the values for the diffusion coefficients, the next step could be to consider a bidimensional and space-variable velocity field, solving the transport equation, and applying other numerical scheme solvers or more sophisticated hydrodynamical models, as one of the limitations that arise when implementing the PDE Toolbox is that it only allows the configuration of an analytical velocity field.
Future studies will focus on calibrating a numerical model by using the optimal coefficients calculated in this study. Velocities and diffusion coefficients calculated at different times will serve as a validation.