During VM generation for the Hikurangi srf provided by Aecom issues were found, meaning that the VM did not contain the SRF fully, even though there are procedures in place in srfinfo2vm to create a boarder that does not contain the SRF.
These issues were primarily related to numerical approximations used that are good enough for most cases, but not for all.
build_corners
The function build_corners(origin, bearing, xlen, ylen) used an iterative method to approach the correct location of the corners. This continued until the distance between the corners was within 10 km of the desired distance.
The new working computes the distance between the origin and the top/bottom mid point and calculates the exact corner from there.
reduce_domain
The function reduce_domain(c1, c2, c3, c4) (c1-c4 are the corners of the VM) would scan from south to north along the edges of the domain and calculate the maximum each side could be reduced by. The origin would then be shifted to match the new center line. The function would then scan from south to north searching for any intersections with land or the srf and bringing in the north and south edges if none were found respectively. The origin was then moved along the north/south line to match the new mid point. During this process the corners were built from the new origin, xlen and ylen multiple times, allowing for the numerical approximations to accumulate into errors.
Issues arose from this method as the distances moved were calculated using the distance at the edge of the boundary, rather than through the middle, which is larger due to the domain being on the surface of a sphere, rather than a plane. The bearing of the mid point was not updated either, resulting in the VM having the incorrect bearing.
The new function takes in the origin, bearing xlen and ylen. It scans from south to north searching for intersections with land or the srf, traveling along the center y line, extending scan lines to the relative east and west. When the amount to reduce is determined the origin is shifted by that amount in a ratio with the mid x length over the edge x length. This results in less error than shifting by the amount without the ratio. The bearing of the new origin is then determined as the right angle from the bearing of the new origin to the old origin
The north/south reduction is then performed, moving the origin along the bearing line by half the difference between the north and south trims, in the direction of the smaller trim. This keeps the origin in the center of the domain. The bearing is then calculated as the direction that is being faced after the origin has moved to account for the new ylen
Results
Using the data from the Hikurangi srf that inspired this review of the velocity model generation code, a map has been created that displays the differences between domain generation between the old and new methods.
The map is available here: https://drive.google.com/open?id=1w95jljHJeafDvT2uA8ZQ9VcvfYUSNCmT&usp=sharing
The old velocity model domain is visible below:
| The srf domain is shown in purple in the center. The initial VM domain is shown in green, and the final VM domain is shown in red. Note that a corner of the final domain is outside the boundary of the initial domain. More importantly a corner of the srf is outside the final domain. This has caused issues causing the computational steps of the workflow to crash, making the simulation useless. |
| The srf is again shown in purple in the center. The initial VM domain is shown in light green, and the final VM domain is shown in orange. Notice that the bottom of the final domain follows the line of the bottom of the initial domain. There is some difference between the two, this is accounted for by the distortion of creating a square on the surface of a sphere defined by the edge lengths. Of primary concern to note is that the srf corners are all contained within the final domain, meaning that this simulation is usable. |
Backing mathematics
For ease of reference, when looking top down at a domain, the edge in the direction of the bearing is the top and the edge in the opposite direction the bottom. The left and right edge follow from this.
In order to accurately determine where the corners of the velocity model domain are we need to know the distance between the center point (origin) and the middle of either of one of the boundary lines.
Without loss of generality for these calculations the distance between the origin and the middle of the top edge is calculated.
For these calculations we use the radius of the earth (noted as R), and assume that the earth is a sphere. Distortion from the earth actually being an oblate spheroid is considered to be trivial.
(1) Domain top down | Here we show a domain as projected onto the surface of a sphere. Note the bulge of the boundary lines to represent that the distance between them is not constant. The value we need in order to calculate the position of the corners is ylen_mid/2, the distance from the origin (red) to the top mid point (light blue). |
(2) Domain profile view | Here is a profile view of the top right quadrant of the domain from the diagram above (Not to scale) with the point at the bottom representing the center of the earth. This is viewed from the right, as the mid right point (purple) is on the left in this diagram. The two values we know are from the purple to yellow point (ylen/2) and from blue to yellow (xlen/2). From this we can get the angular distance θylen/2 = (ylen/2)/R, likewise for θxlen/2. Conversely we only need to determine θylen_mid/2 to determine the length of ylen_mid/2 (red to light blue). By using the radius of the earth and taking the sine of θylen/2we get R*sin(θylen/2), the perpendicular distance between the corner (yellow) and the closest point on the line between the center of the earth and the mid right point(purple). It is useful to note that the plane formed by the four corners, ignoring the surface of the earth, forms a rectangle when bounded by the lines between the corners. Projecting lines from the center of the earth, through the edges of this rectangle and to the surface of the earth give the boundaries of the domain. |
(3) Profile view of the plane following the mid-y line | We should also note two points, the first being the point on the rectangle between the center of the earth and the light blue point, let this be point pymp. The second being the point on the rectangle between the origin and the center of the earth, let this point be point pop. These are both shown on the diagram. The distance R*sin(θylen/2) is the same as the distance between the points pymp and pop. This occurs because the area of the plane bounded by the four corners is a perfect rectangle, and the line projected from the center of the earth through the mid points of the rectangle to the surface of the earth form the mid points of the domain boundaries. We can get the distance between the center of the earth and pympas this is the value R*cos(θxlen/2). The converse for the point pxmp can be observed from diagram (2) as being the distance from the center of the earth to the plane, along the line from the center of the earth to the purple point. Thus, as we have the lengths from the point pymp to the point pop and the center of the earth, the value of θylen_mid/2 is therefore arcsin(R*sin(θylen/2)/R*cos(θxlen/2)). Canceling allows us to remove the Rs from above and below the line giving θylen_mid/2=arcsin(sin(θylen/2)/cos(θxlen/2)). As noted above the the length ylen_mid=R*θylen_mid/2=R*arcsin(sin(θylen/2)/cos(θxlen/2)) = R*arcsin(sin(ylen/2R)/cos(xlen/2R)). Q.E.D. |
In conclusion we have:
