the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 23 Apr 2024
Centroids in second-order conservative remapping schemes on spherical coordinates
Abstract. The transformation of data from one grid system to another is common in climate studies. Among the many schemes used for such transformations is second-order conservative remapping. In particular, a second-order conservative remapping scheme first introduced in 1987 and extended in 1999 to work on the general grids of a sphere has, either directly or indirectly, has served as an important base in a variety of studies.
In this study, the author describes a fundamental problem in the derivation of the method proposed by a pioneer study relating to the treatment of the centroid used as a reference point for the second-order terms in the longitudinal direction. In principle, use of the original formulation may cause damage to the entire remapping result. However, a method's native implementation software includes a preprocessing procedure that tends to minimize or even erase the error as a side effect in many, if not most, typical applications. In this study, three alternative formulations are proposed and tested and are shown to work in a simple application.
- Preprint
(1207 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2024-1101', Moritz Hanke, 24 Apr 2024
Fuyuki Saito found an error in the formulation of the weights for second order conservative remapping paper by P. W. Jones from 1999 (referred to as J99). This error also made it into the SCRIP library, which is based on J99. This software library has been the basis for conservative interpolation in many climate models for many years, which makes the finding presented here a substantive contribution to the community. Unfortunately my limited understanding of the underlying math prevents me from commenting on the correctness of the presented formulas.
The author reproduces the results from J99 and describes in detail the error made therein. This is followed by an analysis of the SCRIP source code and the description of preprocessing step in SCRIP. This step avoids the error having an actual impact on the interpolation results for many common use-cases. Multiple solutions to problem are derived and explained in-depth. Afterwards the impact of error is analyzed in great detail and compared to presented solutions.
Overall the paper is very well written and after a minor revisit, I would recommend it for publishing.
General remarks:
Even through J99 and the SCRIP library have been widely used in the past, the author seems to wrongfully assume that this is still the case for current software (e.g. ESMF). This is apparent in the additional remarks and the summary of the paper. However, ESMF and other software (e.g. YAC or XIOS) use implementations for first and second order conservative remapping, which significantly differ from the SCRIP library and are therefore not prone to the presented error. Second order conservative remapping in these implementations is based on Kritsikis et al., 2017 and not J99.
In the summary of the paper the author encourages the further use of the SCRIP library. However, he fails to mention the various other drawbacks of it, which would lead me to a different conclusion. These drawbacks include inaccuracy for cells close to the poles (also mention by the author) due to how trigonometric functions are used for intersection computation and the misrepresentation of the true cell shapes for everything but RLL rectangular grids (Taylor, 2024). This limits the use of the SCRIP library in my opinion to remapping between two RLL rectangular grids, which could be implemented much simpler and more accurate.
Specific comments:
Line 3: missing reference for "1987" and "1999"
Line 3-4: duplicated "has"
Line 5: missing reference for "pioneer study"
Line 18: Conservative interpolation is just one of multiple methods being used in ESM’s. However, this sentence implies that it is the most commonly used method.
Line 34-37: The CDO’s use code extracted from YAC for the first order conservative interpolation. (see Taylor, 2024)
Line 42: Give Equation number from J99 ("(10)"?).
Section "Introduction": This sections could mention other issues of SCRIP and solutions to this implemented in other software (see General remarks).
Line 73 Eq. (2): "r" is not described
Figure 2: Instead of showing the actual source code, a mathematical description of what the code does might be easier to understand.
Line 245: This may not only happen at the poles. A cell with the longitude bounds of [179;181] may be represented by [179;-179] independent of that latitude.
Line 205: The title could be more concise. In general this paragraph contains in my view too much speculations and opinions of the author. It could probably be shortened without loosing significant information.
Line 256: The implementations of trigonometric function can be more accurate for small absolute values. This may be another reason for this code in SCRIP.
Line 366, 375, 381: repeated use of "naturally...within the cell"
Section "Additional remarks": In this section the discussion of the results for ESMF and other software tested by Valcke, 2022 assumes that it is based on J99, which is not the case.
Line 380: "less impact than a change in magnitude" this is not further explained and quantified. Is this due to numerical inaccuracy of SCRIP for cells close to the poles?
Line 399: "which has the maximum deviation from the pivot longitude within a source cell" this has already been explained above
Line 403: Maybe you could explicitly mention that Figure 4(c) shows deviations from the exact fields, which are independent from the issue discussed in this paper. Additionally, you could state that in order to be able to visualize this issue, you have to compute the difference (d-b/f-b) instead of (d-a/f-a).
Line 522: duplicated "this"
Line 578,585,594: duplication of "https://doi.org"
References:
Kritsikis, E., Aechtner, M., Meurdesoif, Y., and Dubos, T.: Conservative interpolation between general spherical meshes, Geosci. Model Dev., 10, 425–431, https://doi.org/10.5194/gmd-10-425-2017, 2017.
Taylor, K. E.: Truly conserving with conservative remapping methods, Geosci. Model Dev., 17, 415–430, https://doi.org/10.5194/gmd-17-415-2024, 2024.
Valcke S, Piacentini A, Jonville G. Benchmarking Regridding Libraries Used in Earth System Modelling. Mathematical and Computational Applications. 2022; 27(2):31. https://doi.org/10.3390/mca27020031
Citation: https://doi.org/10.5194/egusphere-2024-1101-RC1 -
AC2: 'Reply on RC1', Fuyuki Saito, 19 Jul 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere-2024-1101/egusphere-2024-1101-AC2-supplement.pdf
-
AC2: 'Reply on RC1', Fuyuki Saito, 19 Jul 2024
-
RC2: 'Comment on egusphere-2024-1101', Phil Jones, 21 May 2024
I will first note in this review that my original publication (J99) and subsequent implementation is far from perfect and has some serious issues. I've always wanted to go back and correct those but unfortunately never got the time to do so. I say this to emphasize that the critical review below is not meant as a reactive defense of J99 or the SCRIP implementation. However, in reviewing the paper, I found the author has made some significant errors and incorrect assumptions that negate the conclusion. I do not believe this paper can be published in its present form since it is incorrect.
The first error is the derivation in section 2.1. The author attempts to derive the flux distribution from a Taylor series expansion. However the constraints in equations (4),(5) do not necessarily follow from (3) or at least not uniquely so. They are merely a reasonable and obvious choice among a number of possible solutions. For this reason, neither Dukowicz and Kodis (DK87) nor J99 derive this form from a Taylor series. The two previous papers (JK87, J99) simply show that the flux form:
meets the conservation condition as long as the reference point (r_n) used in the flux approximation is the centroid. It is an assumed distribution that meets the conservation condition. This might seem a minor quibble since the author arrives in a similar place as the two prior papers, but it is important because the author takes the Taylor series approach later as well and this is incorrect.The author correctly notes that equation (6) only holds for Cartesian coordinates. In spherical coordinates, the dot product and the unit vectors are spatially dependent and cannot be formally pulled out of the integral and are more complex in form. The original paper J99 arrives at this form in a different manner. There is a similar issue with the centroid definition, but I will address both of these choices below.
The real problems with the current paper come in section 2.2 where the author incorrectly represents the J99 derivation by stating we use a Cartesian space in lat/lon. This is not the case. All of our derivation occurs in spherical coordinates or a local spherical surface approximation with the appropriate metric factors included. The author can be excused in misunderstanding the derivation since much of the derivation is left to the reader in the original J99 paper. But this mistaken assumption leads to the incorrect form that is the core of the paper.
To elucidate the error, I have to explain how we actually derive the weights in J99 and show some additional steps. We start with the form of the flux approximation shown above (noting again, that this is not a Taylor series expansion):
We also assume the gradient is fixed with the form:
As noted previously the dot product in spherical coordinates is in general not simply a component-wise product as in Cartesian coordinates since the unit vectors can change direction based on position on the sphere. Here we can take two approaches which lead to the same approximation. One is to say that the unit vectors are nearly aligned so that local orthogonality is almost true (it is true for r in this case, but not quite r_n). We use the fact that the local displacement on the unit sphere is:
Then we can approximate the flux as:
The same result can be obtained by using a local quasi-Cartesian approach but including the spherical metric factors. This form is very close to the author's equation (11) except that the first term is the mean flux and the author's pivot point cannot be an arbitrary pivot point. It is required to be the centroid r_n. These differences again arise from the mistaken use of a Taylor expansion rather than the flux form.At this point, the author makes the mistake at the core of the paper. The author assumed we were working in some sort of Cartesian space in (theta,phi) with none of the metric factors and then assumes that a density factor is required in equation (12) to correct the integrals and must occur in all integrals. In fact, we are working in spherical coordinates where the area element dA is defined as
So we do not need this imagined sigma density in equation (12) and equation (18) is incorrect in the denominator.In J99, we compute the centroid, using the standard definition
This leads to the correct latitude centroid in equation (17). The position vector r must be dimensionally (and metrically) consistent with the displacement equation above, so the centroid term in longitude includes the cos(lat) metric scale factor and
Substituting this form of the centroid, we obtain equation (28) (equation 10 in J99). The author's equation 30,31 are incorrect since they depend on the incorrect equation (12). As another aside, we note that a more careful computation of a real centroid should follow, for example, the approach Du, Gunburger and Ju (2003, SIAM J. Sci Computing) in which the centroid is the full 3d centroid constrained to the spherical surface. The centroid here is a very close approximation to that form and is consistent with the dot product assumptions made in the flux approximation so that actual conservation is ensured in practice.The author's error in equation 12 renders the remaining discussion essentially moot. However, I will note that in the later discussion the author misunderstands the longitudinal correction shown in Figure 2. This correction is necessary to avoid computing a line integral across the branch cut in the multiple-valued longitude as explained in section 3(d)(1) of J99.
I hope this review is clear and helps the author understand what was done for J99 and the SCRIP implementation. As noted at the beginning, there are other problems with SCRIP, especially the parametric form for cell sides used in equations 16, 17,18 of J99 and in the SCRIP implementation. Both the choice of a linear form and the failure to include the cos(lat) metric factors in the longitude during intersection computations (even though I did so everywhere else in the paper) are the source of most of the problems in SCRIP and are especially amplified in the 2nd-order terms, making SCRIP difficult to us for higher-order cases.
Citation: https://doi.org/10.5194/egusphere-2024-1101-RC2 -
AC1: 'Reply on RC2 (a brief question to confirm one formulation)', Fuyuki Saito, 21 May 2024
I really appreciate the careful review by Prof. Jones. The derivation you shown helps me understand the detail of the algorithm, which I failed to catch up.
Before my complete response, please let me confirm the formulation of the last equation in your review.
I am still confused with the term on the left hand side: is the latitude (theta) the centroid latitude (i.e., theta_n) ? Or a variable to hold arbitrary latitude (theta, as is)? Or, rather, cos(theta)phi is unbreakable (i.e., the whole terms are for the centroid, as [phi cos theta]_n)?
I suppose only the second one (theta as is) can satisfy Eq.10 in J99 and 5th equation in the review, but it is strange to have a variable latitude term in the centroid longitude. It means that centroid longitude depends on the latitude of a point in a cell to compute the flux approximation.
Another possibility is that the final term in parenthesis in 5th equation to replace with {cos(theta)phi - [cos(theta)phi]_n}, which corresponds to the third possibility above. But in this case the local displacement dr should be dr = d(theta) <theta> + d(phi cos(theta)) <phi> (<theta> and <phi> are the unit vector).
I am afraid that it is a trivial derivation, but I really want to know where I have missed to follow the original scheme.
The other derivation, except for the last, is really clear to me. Thanks a lot for your review. I will come back to this issue in the complete response.
Citation: https://doi.org/10.5194/egusphere-2024-1101-AC1 -
RC3: 'Reply on AC1', Phil Jones, 21 May 2024
Ah, yes. Did get a bit sloppy/inconsistent there. I probably should have stuck with the position vector here (phi only in the centroid) and included the cos(theta) metric only when computing distance as in the flux expansion. Guess I got used to tossing in the metric factor everywhere and it does make for a cleaner weight calculation and is probably better behaved. In the end, I suspect it may not make a large difference - basically the difference between average of the product and product of the average. But better to be consistent.
Citation: https://doi.org/10.5194/egusphere-2024-1101-RC3 -
AC3: 'Reply on RC2 and RC3', Fuyuki Saito, 19 Jul 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere-2024-1101/egusphere-2024-1101-AC3-supplement.pdf
-
AC3: 'Reply on RC2 and RC3', Fuyuki Saito, 19 Jul 2024
-
RC3: 'Reply on AC1', Phil Jones, 21 May 2024
-
AC1: 'Reply on RC2 (a brief question to confirm one formulation)', Fuyuki Saito, 21 May 2024
-
RC4: 'Comment on egusphere-2024-1101', Vijay Mahadevan, 30 May 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere-2024-1101/egusphere-2024-1101-RC4-supplement.pdf
-
AC4: 'Reply on RC4', Fuyuki Saito, 19 Jul 2024
The comment was uploaded in the form of a supplement: https://egusphere.copernicus.org/preprints/2024/egusphere-2024-1101/egusphere-2024-1101-AC4-supplement.pdf
-
AC4: 'Reply on RC4', Fuyuki Saito, 19 Jul 2024
Data sets
Resources of Saito (submitted to GMD) — software and experiment data archives Fuyuki Saito https://doi.org/10.5281/zenodo.10892796
SCRIP-p (p is for pivot) F. Saito https://github.com/saitofuyuki/scrip-p
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 414 | 85 | 32 | 531 | 10 | 10 |
- HTML: 414
- PDF: 85
- XML: 32
- Total: 531
- BibTeX: 10
- EndNote: 10
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1